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1. INTRODUCTION 


To study the effectiveness of various control system design methodologies, the NASA Langley Research 
Center initiated the Benchmark Active Controls Project. In this project, the various methodologies will be 
applied to design a flutter suppression systems for the Benchmark Active Controls Technology (B ACT) Wing 

(also called the PAPA wing). Eventually, the designs will be implemented in hardware and tested on the BACT 
wing in a wind tunnel.. 

This report describes a project at the University of Washington to design a multirate flutter suppression 
system for the BACT wing. The objective of the project was two fold. First, to develop a methodology for 
esigntng robust multirate compensators, and second, to demonstrate the methodology by applying it to the 
design of a multirate flutter suppression system for the BACT wing. 


The contributions of this project are 

1 ) Development of an algorithm for synthesizing robust low order mritirate control laws. The algorithm 

is capable of synthesizing a single compensator which stabilizes both the nominal plant and multiple 
plant perturbations. 

2) Development of a multirate design methodology, and supporting software, for modeling, analyzing 
and synthesizing tnultirate compensators. 

3) Design of a multirate flutter suppression system for NASA’s BACT wing which satisfies the specified 
design criteria 


This report describes each of these contributions in detail. Section 2.0 discusses our design methodology. 
Section 3.0 details the results of our multirate flutter suppression system design for the BACT wing. Finally, 
Section 4.0 presents our conclusions and suggestions for future research. 

The body of the report focuses primarily on the results The associated theoretical background appears in 
the three technical papers that are included as Attachments 1-3. Attachment 4 is a user’s manual for the 
software that is key to our design methodology. 
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2. A METHODOLOGY FOR DESIGNING 
MULTIRATE COMPENSATORS 


2.1. OVERVIEW 
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As we will later see, multirate systems which satisfy our assumptions are periodically time- vary mg. To 
emphasize their periodic nature we will use a double index notation lor the independent variable of a sampled or 
discrete signal. For example, given a continuous signal v ( r). v(m.m represents yU) sampled at the time 
t = tmN + n)T: where the integer N is the period of repetition; T is the sampling period; m =0. i. ... and 

«=0. 1. ... A-l. 

The design methodology presented in the following sections provides tools to model the closed loop system 

in Fig. 2.2. to compute optimum values of A,. B z , C. and D z , and to analyze the performance of the closed-loop 
system. v 

13. MODELING A MULTIRATE SYSTEM 

Two useful modeling tools are the Generalized Multirate Control Law Structure (GMCLS) and the 
Equivalent Time-Invariant System (ETIS). 

13.1. The GMCLS 

The GMCLS is a control law structure which describes a multirate compensator of arbitrary dynamic 
order, with an independent sampling rate for every compensator input, and independent update rates for every 
processor state and compensator output. A multirate compensator with the GMCLS is shown in Fig. 2.2. In 
this figure each element of the continuous plant output v, is sampled at an independent rate. The sampled value 
of v„ y. is combined with the current processor state vector. 5. using the state space structure shown in the 
figure. Each element of the processor state vector, z. is updated at an independent rate. The continuous output 
from the compensator, represented by the vector u. is formed by holding the output from the digital processor, 
with a zero-order-hold. Each element of the vectoi u can be held at an independent rate to form u. 
Conceptually, one can divide the multirate compensator into two parts, the ‘sampling schedule" and the 
digital processor gains. This is the approach used in the GMCLS. The "sampling schedule" is a description of 
when each compensator input is sampled and when each compensator output and processor state is updated, 
while the digital processor gains determine the dynamics of the digital processor. 

2.3.1. 1 Sampling Schedule for a GMCLS 

In general, the sampling and updating of the elements of y,. z. and u in Fig. 2.2 can occur at any time. 
However, to conform to the GMCLS. we require that these sample and update acuvities occur only at integer 
multiples of some fixed time, called the shortest tme period (STP). The actual value of the STP is arbitrary, but 
it is often a function of the hardware and software used to implement the control law. We also require that the 
sampling and updating activities of the sensors, states and outputs repeat themselves after some fixed period of 
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Figure 2.3 Example Sampling Schedule 


Figure 2.4 Aperiodic Sampling Schedule 


time. (This requirement disallows, for example, a system whose sampling period is a function of the time 
require to execute the control software which might vary with control inputs values. ) The penod of repetition of 
the sampling schedule is called the basic time period (BTP). Finally, we define 

Btp 

the integer N = and the value T = STP (2.1) 

In our double index notation, the first index ( m) in, for example, indicates the integer number of BTP’s 

which have elapsed when the sample/update occurred and the second index (n) indicates the integer number of 
STPs which have elapsed within the current BTP when the sample/update occurred. 

Wc can represent the sampling schedule for the multirate compensator graphically, as shown in Fig. 2.3. 
The figure shows a time line for each sampler, processor state, and zero-order-hold. The time line is divided 
into one STP increments. On the left side of the time line is a description of the signal or state being sampled or 
updated. On the right side is a descnption of the particular activity represented by the lime line, e.g., state 
update, sampler, or zero-order- hold. Circles on each time line indicate when a sample or update activity 
associated with that particular signal or state takes place. Usually the sampling schedule is shown lor only one 
BTP since the sampling schedule repeats itself every BTP. 

In most applications, the sampling/updating activities for a given sensor, output or state will be penodic 
within the BTP, as is shown in Fig. 2.3. However, the sampling/updating activities do not have to be penodic 
within the BTP. The only requirement is that the sampling/updating activities have some penod of repetition 
(the BTP) and that they occur at integer multiples of the STP. Once the STP and BTP have been selected, the 
designer can arbitrarily specify sampling/updating acuvitics at any multiple of the STP within one BTP. An 
example of a multirate sampling schedule in which the sampling/updating activities are not penodic within the 
BTP is shown in Fig. 2.4 A sampling policy like this might be used to multiplex multiple inputs through a 
single analog to digital converter. 

2.3.1. 2. Diktat Processor Cains 

The processor gains are the values of the matrices A : , fl : . C z , and D z in Fig. 2.2. Like the sampling 
Nchcdule. they can be periodically time-varying with a t enou of repetition of one BTP. Generally, these 
matrices are free design parameters which can be adjusted by the designer to improve the performance of the 
multiratc compensator. The synthesis algorithm discussed in Section 2.4 erm be used to calculate optimum 
values lor these gains. 
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2. 3.1.3. State Space Formulation of the GMCLS 

A compensator with the GMCLS can be modeled as a periodically t.me-varying discrete-time svstem. The 
state space form o t the GMCLS is given by 


c(m,/i+l) = Ag(n)z(m,nl + B g (n)\imjt) (2.2a) 

u(mj i)= Cg(n)zinin)+ Dg(n)z(m,n) (2.2b) 

where 

z(nun) = U(w,/z) T x(mM) J u(m,n) T \ r (2.3) 

and is used to model the sample and hold activity from u(m.n) to u< m.n). The form of A„. C and D e 

is given in (Berg. Mason & Yang 1991 ] and (Mason & Berg 19921 which are included as Attachments I and 2. 

We should emphasize that Eqn. (2.2) is used to model the complete sampling/updating activities and 
dynamics ot a multirate compensator. It would not be used in the actual implementation of the compensator. 
When implemented, the sample and hold activities of the inputs and outputs would be performed by appropriate 
hardware. The only dynamics to be calculated are those associated with the processor state vector z. 

2. 3. 1.4. Factored Form of the GMCLS 

Equation (2.2) is a convenient form to model the general multirate compensator. The difficulty with 
Eqn. (2.2) is that it ties up the digital processor matrices, A z (n), B : (n). C-(zi), and D z (n), in the model matrices 
.4 ? («), Bgin), C g {n), and D g in). The matrices A z (n), B z (n), C z (n). and D.in). which describe the dynamics of 
the digital processor, are the unknown design parameters which we will later optimize. We can separate the 
processor dynamics matrices from the model matrices as follows. 

Define the composite compensator matrix: 


and factor Pin) as follows 



L Bgin) 


Cgin ) 
Agin) . 


(2.4) 


Pin)= .Si (n)P z (n) S 2 (n) + S 3 (n) 


(2.5) 


where P z (n) = 




D z (n) 


B-in) 


C z (n) 
A z (n) . 


(2.6) 


and S i , 5 2 and S 3 are the switching matrices defined by the sampling schedule for the compensator. Their 
exact form is given in (Mason 1992] and (Mason & Berg 1992] (see Attachment I ). 

It is important to note the difference between Pin) and P z (n) in Eqn. (2.5). Pin) is a periodically time- 
varying matrix defined by Eqn. (2.4) It includes all the information about the processor gains and the 

sampling/update schedule. P z (n) contains only the gains for the processor dynamics and is independent of the 
sampling schedule. 

2.3.1. 5. Implemer, tation 

The Generalized Multirate Control Law Structure (GMCLS) provides a framework for dealing with 
multiple sample/update rates, time delays, and periodically t.me-va' ing gams in a digital control system. It 
gives the designer freedom to either sd t the "sampling schedule" that best solves tU prot'icin, or if necessary, 
to use the sampling schedule" dictate^ listing hardware and software, with out lias mg to worry about the 
bookkeeping involved with multiple rates i time delays. 
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Ir practice, the GMCLS is implemented in software and is rarely used directly by the designer. The 
designer need only supply the sampling schedule and values for the digital processor gains to provide a 
complete compensator description. This description can then be transformed directly into a single-rate 
periodically time-varying system using the GMCLS. 

The GMCLS is used extensively by the synthesis algorithm described in Section 2.4. and by the modeling 
and analysis software referred to in Section 2.5. Documentation for this software is provided in Attachment 4 

2J.2. The Equivalent Time-Invariant System (ETIS) 

A multirate compensator with the periodically time-varying structure discussed in Section 2.3. 1.3 can be 
further transformed into a single-rate Equivalent Time-Invariant System ( ETIS) with the form shown below 

.tfm+I.O) = A E xim, 0) + Bpugm, 0) (2.7a) 

\£<»n,0) = C[rMm,0) + DgUpim, 0) (2.7b) 

where 



Vy(m,0) 


u(m,0) 

ygi tn ,0 ) = 

Vj(m.l) 

und i/£<m.0) = 

u(mA) 


_v 5 (mjV-l \ m 


-w(mjV-l)- 


We use the subscript £ to denote vectors and matrices strictly associated with the ETIS. See |Meyer & Burrus 
1975] or (Mason 1992] for a definition of Af, fl/r, Cg and D f. 

A key feature of an ETIS is that a multirate, or periodically time-varying system will be staole if and only if 
its ETIS is stable (Kono 1971], Also notice that the ETIS input/output vectors are composite vectors 
containing the mput/output values of the multirate (or penodically-time varying) system at N sampling times. 
Consequently, an ETIS is always MIMO even if the original system is SISO. If the multirate system has p 
inputs, q outputs and a sampling period of one STP then the ETIS is a single-rate linear time-invariant system 
with N p inputs, N q outputs and a sampling period of one BTP. 

2. 3.2.1. Implementation 

"^ ie f TIS is fundamental to the analysis of multirate systems. It allows one to evaluate the performance 
and stability of complex systems comprised of multirate, periodically time-varying and/or single-rate 
components using only techniques developed for linear time-invariant single-rate systems. For example, to 
evaluate the stability of the system in Fig. 2.2. we would first transform the multirate compensator into its ETIS 
with a given value for N. Then we would discretize the plant at the STP of the compensator using a zero-order- 
hold and transform the resulting single-rate system into an ETIS using the BTP of the compensator. Next, the 
plant and compensator ETIS s could be combined in feedback just as if they were traditional single-rate 
systems. Finally, we could determine the stability of the original multirate sampled-data system from the 
eigenvalues of its closed-loop ETIS. 

Documentation for software capable of transforming multirate and single-rate systems into their ETIS’s is 
provided in Attachment 4. 

2.4. SYNTHESIZING A M ULTIRATE COMPENSATOR 

When designing a multirate compensator for the system in Fig. 2.2 there are three components one must 
consider: the compensator structure (this includes the dynamical order of the digital processor), the sampling 
schedule, and the values lor the digital processor gains. In our design methodology the compensator structure 
and sampling schedule are selected by the designer based on the open- loop plant dynamics, the hardware 
constraints, if any. and the desired closed-loop performance. Values for the digital processor gains are 


s 


calculated by our synthesis algorithm so as to provide optimum closed-loop performance for the chosen 
compensator structure and sampling schedule. In the following paragraphs we discuss compensator structure 
and sampling schedule selection, and provide a brief description cl our synthesis algorithm. A complete 
discussion of the algorithm is provided in Attachment 4 and in [Mason & Berg 1992] (also included as 
Attachment I ). 

2.4. 1 . Compensator Structure and Sampling Schedule Selection 

The choice of compensator structure and sampling schedule is problem dependent. It depends on the 
hardware constraints, the open-loop plant dynamics, and the design objectives. Two often used multirate 


compensator structures are worthy of mention, however. They are successive loop closure and coupled 
successive loop closures. (Also see [Berg 1986) for a discussion of successive loop closures.) 

2.4. 1. J. Successive Loop Closures Structure 

The simplest multirate compensator structure is successive loop closures (SLC). This structure consists of 
multiple decoupled single-rate control loops, each loop operating at a unique sample/update rate. The state 


space representation of a SLC structure with two loops is 

1 y slow t n J j 


x slow ( n + I 1 J 


a fast 

0 


Jr 


m) 


a slow J ( x slow ( ** ) 


"fast 'I 
0 b slow 


(2.9a) 


“fasting Jc fast 
U stow( n )) 0 


0 

c slow , 1 x slow ( ** ) J 


dfast 

0 


0 

dsiow 


>' fasti"*) 1 

ysiow ( n ) J 


(2.9b) 


where v represents the sampled input from the sensor and u is the output to the zero-order-hold. The subscnpts 
fast and slow denote inputs, outputs and states which are sampled/updated at a fast or slow rate, respectively. 

SLC is best applied to control problems where the closed-loop dynamics are comprised of some fast and 
some slow dynamics with the bandwidths of the two separated by at least a factor of four. In this type of 
problem, the “fast” loop(s) of the SLC compensator, operating at a fast sampling/update rate, would be used to 
control the high bandwidth dynamics, while the “slow” loop(s), operating at a slower sampling/update rate. 

would be used to control the low bandwidth dynamics. Problems such as these usually fall into one of two 
categories. 

In the first, the open-loop system exhibits both fast and slow dynamics. The multirate compensator is used 
to improve the performance of this system without drastically changing the fast or slow bandwidths. An 
example of this type of problem is an aircraft yaw damper/modal suppression system. The aircraft is open-loop 
stable and has some fast dynamics associated with the flexibility of the airframe and some slower dynamics 
associated with the yawing motion of the entire aircraft. A multirate compensator for such a system might 
consist of a high bandwidth loop to damp the airframe vibrations and a low bandwidth loop to improve yaw 
damping. 

In the second type, the open-loop dynamics of the plant are arbitrary, but in feedback with the compensator 
the closed-loop system exhibits the characteristic fast and slow dynamics. These systems usually have a 
decoupled structure where sets ol open-loop modes are strongly controllable and observable with a particular set 
ol inputs and outputs and weakly controllable and observable with the remaining inputs and outputs. An 
example of thi. type ol system is the two link robot arm (TLA) used in [Berg, Amu & Powell 19881, and in 
[Yang 1988). All four of the open-loop poles of the TLA are at the origin of the “s" plane. The plant has two 
inputs and two outputs. Only two of the modes can be controlled with any one input. Similarly, only two of 
these modes can be observed with any one output. In the multirate design, one input/output pair is used to place 


9 


two ot the closed-loop poles at a high frequency and the other input/output pair is used to place the other two 
closed-loop poles at a low frequency. 

Sample rate selection for the individual control loops of a SLC design follows the same guide lines used in 
single-rate sample rate compensator design: the sample rate for each SLC loop should be 5 to 20 times faster 
than the closed-loop bandwidth desired for that loop. See (Franklin Powell & Workman 1990] for a discussion 
ot sample rate selection tor single-rate systems. 

2.4. 1. 2. Coupled Successive Loop Closures Structure 

The coupled SLC structure is the same as the traditional SLC structure except the designer can include 
cross feed terms which coupie the fast and slow inputs and outputs of the design. In the state space formulation, 
cross coupling is represented by non-zero off diagonal terms in the compensator gain matrices. An example of 
a compensator structure with cros . feed from the slow sampled sensor to the fast sampled/updated control loop 
is given in Eqn. (2.10). 

Klow(« + l)j 

|«/a«(ffi)] _ 

1 U slow( n ) j” 

This structure is best applied to systems which have coupling between their fast and slow closed-loop dynamics. 
See [Yang 1988] for a discussion of cross feed in the TLA problem. 

2.4.2. Optimizing the Digital Processor Gains 

Having chosen an appropriate compensator structure and sampling schedule, the designer can use our 
synthes>‘ algorithm to calculate optimum values for the digital processor gains A z . B z , C z and D z such that the 
closed-lc,,, system in Fig. 2.2 minimizes a quadratic cost function. 

The primary design parameter for the synthesis algorithm is the quadratic cost function. By selecting an 
appropriate cost function, the designer can influence the performance of the resulting closed-loop system. The 
cost function minimized by our synthesis algorithm has the fori.; 


11 fast 0 
^ a slow 
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(2.11) 


where J is the cost associated with the closed-loop system shown in Fig. 2.2. The vector y c is the continuous 
criterion output and u is the continuous control input. (?i , Qi and M are the cost function weighting matrices 
and are free design parameters. 

The cost function in Eqn. (2.11) has the same form in a continuous time LQR design. Thus the cost 
associated with the optimized multirate compensator and that of an LQR design can be compared directly. The 
designer can also use this fact to help select appropriate values for Q \ , Q2 and M. 

To improve the robustness of the compensator, the synthesis algorithm can optimize the digital processor 
gains for multiple plant conditions simultaneously. The resulting compensator will stabilize the each plant 
condition and provide overall optimum performance. This is accomplished by minimizing the new cost 
function of Eqn. (2.12) which is the sum of the costs associated with each plant condition. 



H ![><■,( oil 

H , . .«,<0_ f 


( 2 . 12 ) 
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Here 7, is the cost associated with the i"« plant perturbation and there are Np plant perturbation* 
Optimum values of A z , B z . C z , and D z . occur when 

_0, ^5T = 0 - = () ’ and ~ 0 or equivalently when -£y- = 0 

‘ ' Or-, 


( 2 . 13 ) 


3) 


Eqn. (2 g ?3Mo det* * typC numcncal search and a closed form expression for the gradients in 

satisfied R f e ^,' ne ^ ^ ^ ^ d ' gUa * P rocessor 8 ains such that the conditions in Eqn. (2.13) are 

Eqn. 2J3 T^e Iv!^ 0 " V*"* ^ ***"*« 1 for a closed fo ™ expression for the gradients in 

processor gains and th ^ *" ‘ teraUVe prOCess t0 determ »ne optimum values for the digital 

processor gains and the user must prov.de the software with an initial guess for A.. B z , C-. and D-. The initial 

guess must stabilize every plant condition considered in Eqn. (2. 12). 

2.4.3. Implementation 

In practice, the steps tor designing a compensator with our methodology are 

1 ) Construct a continuous LQ regulator for each plant condition which achieves the desired performance 
tor that condition. 

21 B ““ d 00 lhe des,red ck »««°°P dynamics and the consiramts imposed b, the system hardware, 
choose an appropriate compensator structure and sampling schedule. 

Usmg the chosen sampling schedule and compensator structure, design a compensator which 
stabilizes all plant perturbations. When the desired compensator structure ,s one of the two structures 

“J" ' he P 7 Vi ° US SeCt, ° n ’ thc designer ca " «• successive loop closures to find a stabilizing 
the digital processor gains. In successive loop closures, the plant is stabilized by closing one 
oop at a time trom one set of inputs to one se, of outputs. To obtain a mult, rate compensator each 

°h P H ,S . k USmg 3 d,ffereM Sampling/U f )date rate - Wh en. due either to a complex sampling 
sc e u e. or t e complexities ot the control problem, successive loop closures cannot be used to find a 

s a i utng value tor the digital processor gains, use Yang’s algorithm (see (Yang 1988]). This may 

^"1 rrT 0d T ,Ve at flrSt ’ S,nCC ° ne 0f the reasons for developing our algorithm was the 
computational inetficences of Yang’s algorithm. However, our experience has shown that. ,n 

genera . ang s a gonthm converges to a stabilizing compensator fairly rapidly. It is the computation 
time associated with optimization of this stabilizing solution that tends to be excessive. 

^fi^Tr" 1 Va r UeS f ° r ,hC dlg,,a ‘ Pr ° CeSSOr gamS US,ng ,he s y nthes,s a,gor ' thm of 

h u 1 C ° St " We ' gh,ing matr,Ces ,or the °P ,,ni ‘ z ation are the same as those used to 

design the LQ regulators in Step I. The starting point for the optimization ,s the stabilizing 
compensator designed in Step 3. 

See Attachment 4 fo, the complete docutnenlahon of the softwatc that tmplements the Synth™, algonthm 
2.5. ANALYZING A MULTIRATE SYSTEM 

tradutn'Tr TT " d ‘ ffiCUh ^ ‘ he peri « dic " ature of a * -stem implies tha, a 

til I* 5 "°' CX,St Th *- COmm ° n 3nalySiS ,00 ' S SUch as *«*** res P° nsc or Nyquist 

tZTuZ " T l ° mUU,rate SyS,CmS ° Ur S0luti0n “ transform the mu, tie system ,nto 

S T : Tate SyS,Cm ’ ,hC ETIS - a " d ‘ hen appl > ««*«»■* single-rate analysis techniques 

G X Th 7?°"" ° 1 ' ETIS <NO,C: WC WntC thC Z ' Transform of an ETIS where N=BTP/STP as 

I ■ e 0 0W ] ng para S raphs discuss five useful tools for analyzing the performance and stability of a 

multirate system based on its ETIS. 7 


4 ) 



Figure 2.5. Plant/compensaior configuration Figure 2.6. Piant/compcnsacor with uncertainty 


2-5*1’ Gain and Phase Margins 

In Section 2.3.2 we noted that a multirate system will be stable if and only if its ETIS is stable. Therefore, 
we can determine whether the multirate system is stable by applying the Nyquist criterion to its ETIS. Since all 
but tnvial ETIS s are MIMO. we must use the multiloop Nyquist stability criterion. The multivariable Nyquist 
is a plot ot the eigenvalues of the ETIS loop transfer function as the discrete variable z traverses the unit circle 
[MacFarlane 1 970] [Maciejowske 199C 1 

When the multirate system is SISO we can obtain traditional gain and phase margins from the multiloop 

Nyquist plot. Let G E (z N ) be the ETIS loop transfer function and let A be some constant gain and phase 
uncertainty at the plant input. If 

Mz)= kJ® where k is a scalar (2.14) 

then A (2.15) 
where / is and Nx N identity matrix 

Now the new loop transfer function with the gain and phase uncertainty of Eqn. (2. 15) can be written as 

He (z' V >loop = Gf ) k ( 2 . 16) 

The multiloop Nyquist plot ot //£(c‘ V )| oop is just the multiloop Nyquist plot of G E (j ,V ) scaled by the gain k and 
rotated by the phase shift 6 - the same as in traditional SISO Nyquist plots. Gain and phase margins for the 
multirate system can therefore be obtained from the multiloop Nyquist plot of G E <- V ) bv determining the 

values ot k and 0 which destabilize the ETIS. (See (Thompson 19861 for an alternate derivation using Kranc 
operators.) 

When the multirate system is MIMO, the gain and phase margins calculated by this procedure apply 
simultaneously to all inputs and outputs, and are consequently not realistic measures of robustness. To obtain 

realistic measures ot robustness tor a MIMO multirate system, a norm based approach such as singular value 
analysis is required. 

15.2. Singular Values 

Singular values are useful for measuring the robustness of MIMO multirate ./stems. The key step in 
multirate singular value analysts is transforming the multirate system in Fig. 2.5 into an ETIS system which has 
the output feedback torm shown in Fig. 2.6. Since the multirate system will be stable if and only if its ETIS is 
stable, the closed-loop system in Fig. 2.5 will be stable for a given value of A provided the closed-loop system 
in Fig. 2.6 is stable tor a corresponding value of Af. Thus we can use single-rate techniques to evaluate the 
robustness ot the ETIS system and relate those results directly to the associated multirate system. 







2. 5.2. 1. Unstructured Singular Value Analysis 

A bound on the smallest value of a(A£) for which A£ destabilizes the system shown in Fig. 2.6 can be 
calculated using unstructured singular value analysis. This system will be stable for all A£ such that 

a(A£ ( z N)) < ~ for all z on the unit circle (2.17) 

g(Ge(z n )) 


(see (Maciejowski 1989]). This result, however, is only a measure of the size of the smallest destabilizing A £ 
and is generally not a measure of the size of the smallest destabilizing uncertainty A. Because the input/output 
vectors ot an ETIS are composite vectors, containing the input/output values of the multirate system at N sample 
times, Af can be a complex function of the values of A at N sample times. (The relation between A£ and A is 
given by Eqn. 2.7.) The size ot the smallest destabilizing A£ found using unstructured singular value analysis is 
only a conservative estimate ot the size of the smallest destabilizing A. This estimate accounts for not only the 
fictitious perturbations normally associated with unstructured singular values, but also for time-varying and 
non-causal perturbations. 

Consider the simple case where A is a constant. From Fig. 2.5 we have that 


vi' = A V' 


For an ETIS with N=2 


w E = A e v e or 


w(m,0)l 
w(m,l) J 


A, | Aj 2 fv(m,0)l 
A 2 i A 22 ^ \ v(m,l) j 


(2.18) 


(2.19) 


A destabilizing A£ determined by singular value analysis might, for example, include block diagonal elements 
in A £ which are unequal, e.g. A\ \ * A 22 - This corresponds to a time-varying perturbation because the gain 
between w and v varies with time. Another such A£ could include non-zero upper block diagonal elements in 
A£, e g. A 12 * 0. This corresponds to a non-causal perturbation because a future input, v(m, I), can affect the 
current output w(m,0). 

We can eliminate this conservativeness by restricting the allowable perturbations in A £. This leads directly 
to structured singular value analysis. 


2* 5. 2. 2. Structured Singular Value Analysis 

In order for the ETIS uncertainty A E i represent the actual uncertainty A, its structure must obey 
Eqn. (2.7). Finding the size of the smallest destabilizing A £ subject to Eqn. (2.7) requires the solution of a 
structured singular value problem. For the system in Fig. 2.6 we define the structured singular value, ji, as 


mo E (z N )) = 


| 0 if det(/ - G E (z N ))* 0 for all A e Ag£> 

min (ct(A(j)) 1 such that det(/ - C E U‘ V )A £ u‘ V )) = 0 
AeA fl0 


-1 


otherwise 


( 2 . 20 ) 


where Abd < s the form of the permissible block diagonal perturbations A and the structure of A E must satisfy 
Eqn. (2.7). The size of the smallest destabilizing perturbation CT(A m j n ) satisfies 

■=77 = sup n(G E (z '' )) where ; v = e j0 (2.21) 

a(A min ) « 


For a discussion of ft and A BD see (Doyle 19821- 
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A[,(z“)- j 


A(c) +- A(-z) 

z(Mz)-M-z)) 


.-1 


(AiZ)- A(-z» 
A(z) + A(-z) 


( 2 . 22 ) 


to 2 "' 7 7* !0lVe Eqn ' l2 ' 201 Wl ' h <0 ta« 

.1* uncertainty, a is , cons* aT s ,£l 7 ***•" "**»<* *■•*» When, however, 

wi* a repeal block “ ETIS “— »• 4* » - • eons* 

A £ = dia *<A. ^ A) with N blocks. (2 23) 

method for estimafinr^wh^rA ^ ^ *** M * W0Ck diag ° nal Struc£ure 0ne simple 

<«ing n when A is strictly diagonal is den ved in [Safonov 1 982J. It is 

MG n (z N )) < .nf(mabs(DG v (^) D ->)) = A p (G N (z N )) (2 . 24) 

where abs( > 4 ) is a matrix such that fabsM )) 14 i- *1 . i rh 1 

eigenvalue. J ^ 15 e 1 J e emenl °* A\ a °d is the Perron -Frobem us 

2. 5. 4. 3 Implementation 

The procedure for performing singular value analysis via the ETIS is as follows 
1 ) Transform the problem into the form shown in Fig. 2.5 

.“nlr”' “ ^ ST? ° f “ d «*- - ™ - PI- using the ,Vof 

Combine the ETO of ,he plant tmtl compensate, to obtain the closed-loop system show, Fig. 2.6 

de7b7i^uKe^n,y'l7' S, " 8Ular '" l “' b “' d ana ' 1 ' sls 1001 10 com P u " «“ »« of the smallest 

"1'““ ET,S mcenamy i£ 

might be conservative unless shoe, wed ^ - ' * — 

Maximum RMS Gain 

as a^d~ ™?,r, zz:Tr s>s " m ,s ,he - - - *—•. ■* ^ 

RMS gatn of a SISO m.“ !y£m IT 7 “ "“"™" 5 ' SKm H °— *• martmum 

value, shown in Eqn. 12 25) plays the same ml 7' " “ °' ETIS l,amf " f »"otion. This 

q (*..*.5). p |a y s ‘he same role as the maximum Bode plot gam of a single -rote system. 


2 ) 

3) 

4) 

5) 


2.5.3. 


sup ^ S(v(ww> > = _ RMS(v f (m.0)> 

RMS<«)#o RMS(u(m.n)) rms«£ ,» 0 RMS(«f(m.O)) “ M Ge(z )H * 


(2.25) 


rms ga,n ,° f s,so ° r m,m ° — * 
work oflStvashankardt KhLglekt iW, * “ * "* ° f ° £ *" "» "* »l— 

muu“^:r,~7::,r ,e 7' r ^ *- .» - 

' °* s no ' necessarily have the simple form stnltu Tm). Instead,, „ comprtwd of the tum of 
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sinusoids o t several distinct frequencies. Details on computing the signal of maximum RMS gain for a 
multirate system are given in [Mason & Berg 1992] (Attachment 1). 

2. 5. 3. 1 Implementation 

One simple method for determining the Hoc norm is to plot the maximum singular value of Gg as z 
traverses the unit circle. H UGe) is then the peak value on that plot. 

It is important to remember that Eqn. (2.25) is a measure of the discrete RMS gain between the discrete 
inputs and outputs of interest. Often the designer is interested in calculating the maximum RMS gain between a 
continuous input and output of a sampled-data system. A good estimate of the RMS gain in this case can be 
found by sampling the continuous input and output of interest at a fast rate. The result is a multirate system - 
the input and output of interest are sampled/updated at a fast rate while the other inputs and outputs are sampled 
at the rate appropriate for connection to the multirate compensator. (This is also useful for determining the 
inter-sample behavior of a sampled-data system.) The maximum RMS gain can then be calculated using the 
ETIS of this new system. 

2^.4. Steady-State Covariance 

A common measure of performance is the steady-state covariance of select outputs in response to a 
disturbance input. In a multirate system the “steady-state” covariance values are periodically time-varying. 

Fortunately, the periodic “steady-state” covariance values at each sample/update time are straightforward to 
calculate using the ETIS. 

It is easy to show that 

£{y(/n,0)>(m,0) T } 

£{y E y£} = £{y(m,l)v(/n,0) r } 

_£{>’( w, jV - 1 )y(m, O) 7 ^ \ 

£{y(m,0)v(/n, 1)^} ••• £{y(m,0)v(ftt, N - 1)^} 

E { y( m, 1 ) v( m, 1 ) T ) E{ y( m , 1 ) y( m, N - 1 } 

£{_v(m, N - l)v(m, 1)^} • £{v(m, N - l)y(m, N - 1)^} 

The diagonal block elements of Eqn. (2.26) contain the steady-state covariance values at each sample/update 
time ol the corresponding multirate system. Therefore, the steady-state covariance values can be found by 
calculating the ETIS of the multirate system and computing the steady-state covariance values of the ETIS using 
the discrete Lyapunov equation. Refer to [Kv/akemaak & Sivan 1972]. Algorithms for calculating discrete 
covariance values are widely available (e.g., in Matlab and in Matrix x X 

2*5.5. Time Domain Simulations 

Time domain simulations are straightforward to compute using the ETIS and Eqn. (2.7). As noted in 
Section 2.5.3, inter-sample behavior can be obtained by sampling the continuous inputs and outputs at an 
arbitrarily fast rate. Documentation for the M-Filc mrsim, which generates a time domain simulation of a 
multirate sampled-data system using the ETIS is provided in Attachment 4. 

2.6. SUMMARY 

The tools presented in this section form the foundation of our multirate design methodology, and provide a 
unified approach to multirate modeling, synthesis and analysis. Using these tools one can model a complex 
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3. APPLICATION OF THE MULTIRATE DESIGN METHODOLOGY TO THE 
DESIGN OF A FLUTTER SUPPRESSION SYSTEM FOR THE BACT WING 

3.1. INTRODUCTION 

To demonstrate some of the advantages of multirate control and the capabilities of our design methodology, 
we designed several flutter suppression systems for NASA’s BACT wing using the methodology in Section 2. 
A summary of our designs is presented in the following paragraphs. In Section 3.2 we describe the model wing 
and its open-loop characteristics. In Section 3.3 we discuss our design goals and constraints. In Section 3.4 we 
discuss our design approach and the details of the design process. In Section 3.5 we present our flutter 
suppression system design results. Finally, we end the chapter vith some concluding remarks in Section 3.6. 

3.2. THE MODEL WING AND ITS OPEN -LOOP DYNAMICS 

3*2.1. Model Wing Description 

The BACT wing was developed by NASA Langley for the Benchmark Models Program. It consists of a 

rigid airfoil mounted on a flexible base. The base, called the Pitch and Plunge Apparatus (PAPA), provides the 
two degrees of freedom needed to model classical wing flutter. Our designs used the single control surface (CS) 
located on the trailing edge of the airfoil and two accelerometers, one near the trailing edge (TE) of the airfoil 
and one near the leading edge (LE). A diagram of the BACT wing is shown in Fig. 3. 1 . A detailed description 
of the BACT wing can be found in [Durham, Keller, Bennett & Wieseman 1991] and [Bennett. Eckstrom, 
Rivera. Dansbeny, Farmer & Durham 1991], 

The flutter suppression system was designed using a 16 ,A order linear state model of the BACT wing 
developed by NASA Langley’s Structural Dynamics Division. This model consists of 4 rigid body states 
corresponding to the pitch and plunge modes, 6 unsteady aerodynamic states, a second order actuator model, a 

second order Dryden filter, and two first order anti-aliasing filters. A block diagram of the mathematical model 
is shown in Fig. 3.2 on the following page. 

We were provided with 24 different mathematical models of the wing. These models describe the motion 
of the wing in freon at 24 different operating points. The operating points include dynamic pressures above and 
below the cnucal flutter pressure at three different mach numbers. See Table 3. 1 on the following page for a 
summary of the operating points. 



Figure 3.1, BACT wing 
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Figure 3.2. Block diagram of B ACT wing 


Table 3.1. Operating points for B ACT wing. All operating points assume Freon medium 


Dynamic Pressure <psf) 

(Nominally unstable operating points are in gray) 


Mach 0.50 

75 

100 

122 

132 


175 


225 

Mach 0.70 

75 

100 

125 

136 

146 

175 

20O 

225 

Mach 0.78 

75 

100 

125 

141 

151 

175 

2ti 

225 


3.2.2. Open- Loop Dynamics 

The response of the open-loop B ACT wing model at each operating point is characterized by two dominant 
modes - the pitch and plunge modes. The poles associated with pitch and plunge at mach 0.5 and 75 psf are 
indicated on Figs 3.3-3.4. As the dynamic pressure increases, one pair of these dominant poles moves towards 
the right half plane and eventually crosses the imaginary axis at the flutter stability boundary. Figures 3.5-3.6 
show the migration of these dominant modes as dynamic pressure increases. The locations of the open-loop 
poles not shown in the figures remain relatively constant. 

The dominant pitch and plunge modes are observable at all operating points with either the TE or the LE 
accelerometer outputs and are controllable at all operaung points using the CS command input. The zeros of the 
CS command to TE accelerometer and the CS command to LE accelerometer transfer functions are shown in 
Figs. 3.3-34 for an operating point ol mach 0.5 and 75 psf. As dynamic pressure increases, the non-minimum 
phase zeros associated with the TE accelerometer migrate into the left half plane. The minimum phase zeros 

that are associated with the LE accelerometer and located near the dominant poles migrate into the right half 
plane. See Figures 3. 5-3.6. 


At low dynamic pressures the transfer functions from CS command input to both the TE and LE 
accelerometer outputs are non-minimum phase. Non-minimum phase systems are typically more difficult to 
control than minimum phase systems. An alternative output is one which measures the difference between the 
two accelerometers. This new output is essentially pitch acceleration. The CS command to pitch accelerauon 
transfer lunctton is minimum phase for all operating points. Figure 3.7 shows the locations of the zeros near the 
pitch and plunge modes as dynamic pressure increases. It turns out that the BACT wing is fairly easy to control 
using this new output. The problem is that the pitch acceleration output is artificially created and assumes 
perfect measurement of TE and LE accelerations. In reality there is some uncertainty in the TE and LE 
acceleration measurements that must be accounted for in any design. Therefore we did not use the pitch 
acceleration output directly in our designs. We did. however, use the pitch acceleration output to determine an 
initial stabilizing compensator for the synthesis algorithm. This is discussed further in Section 3.4.3 














3J. DESIGN GOALS AND CONSTRAINTS 

wjn 4 ° f ,hC dCSign WaS '° SymheS,ZC 3 mUh,raU fluMcr su PP ress '°" system which stabtlizes the BACT 
g operating points. In addition to stability. NASA Langley specified the following constraints. 

Control Activity Constraint: For unity RMS white noise input disturbance < 1 in/sec RMS), the steady- 
. ‘ covariance of the CS deflection must no. exceed 0.0625 deg^ (0.25 deg RMS), and the CS 

deflection rate must not exceed 65 deg^/sec^ (8.0 deg/sec RMS). 

* 3 ',7 Rate * CStna '™ s 1716 mimmum sampling period is 0.005 seconds. Fo, multirate sampling all 
sampling periods must be multiples of 0.005 sec. P * 

Compmmioml D,hy : All compel,*,,,,, mm, * d« S ,g ne d wrth a minimum U.005 »d comp u ,a„„„,| 


I 



21 


and IT thC COm P cnsator *»P*‘ and output must be ±6db 

margins based on the singular value The sne r h° Sem ° TS ' WC US€ ,he general ' 2ed gain and phase 

— or 0.75 f „ r *. ma c;:;: 8Ul ^r:fi 8 r r pha ” n,a,8,ns » • ■— 

(see [Mukhopdhyay & Newsom 1984]) P ,CaUVe uncer *ainty at the compensator inputs 

34. FLUTTER SUPPRESSION SYSTEM DESIGN 

stops for this design were d ' SCIIS ^ d S ““°" 2 10 desl 8" Holler suppression syslem. The specific 

' ‘ SSTSS S “ Ch *“ *• BACT Wi " g ta “ - - tQ regulator saUsfws toe 
£jT an a Ppro p,i M e nmhirato compensator stoic, urn and sampling schednic based on tois LQR 
Find a se. of processor gains so ton, die compensator stabilizes toe BACT wing 

Clieck toe perfoimance and robustness of toe dosed-loop system 
Iterate on items ( 1 )-(5) as required 
We elaborate on the details of each step in the following paragraphs. 

3.4. 1. Selecting the Cost Function Weights 

minimizing a quadratic ^ ° f ' he c ° rn P ensa,or ’s digital processor gains by 

mult, pie plant conditions simultaneously We^sed ih , T optimizatlon *>e performed for 
ensure that the compensator stables ° f «*■»• to he' < 

points for the optimization we used six representative °P* rat,n g P° lnts Instead of using all 24 operating 
of mach number and dynam.c pre^Td Two! i SU ^ ^ 0pe ™" g **« « ■*« extreme! 

operating points are listed in ^ betWee " the extremes. These 

Section 3.4.2 4 we included four additional operation TT P38e ‘ J* ^ faU “ t0lCraiU dCS ' gn discussed in 
Table 3.2. 0pCratmg P °' n,S a ‘ mach 05 ° These operating points are grayed in 

•Hie weig^ts^were'based ^ona! c^ntinmnjs^LQlTd^^ fights for the synthesis algorithm s cos, function, 
command input of the BACT wing. The cos, function h^'the fom ^ ** P “ Ch P,U " ge ^ ‘ he CS 


2 ) 


3) 

4) 


5) 

6) 


^ = f lim£{x(r) T Gl , (r) + tt( ,) T Q2l4(/ )} 


(3.1) 

modalized version of the BACT wing modeT States jr” aSS ° C ' a,cd wi,h the P' ,ch and P |u "ge mode in a 
migrate to the left as dynam.c pressure mcreases. se'eFig^ 77^° "*. C ° mpleX Con ^ M **• 
conjugate poles which migrate to the risht as dv„» mi z. 3 and ** corres P° nd »° the complex 

cause instability in the BACT wing a igh dtnaTc “ mcreases. see Fig. 3.5. The latter set of poles 

For each operating pomt. the n 7Zo ^ ^ CS C ° mmand sig " al 

and plunge modes was greater than 0.07 and the RMS^onu I **" S<> ^ *** closed ' loo P doping of the pitch 

For comparison, the damp.ng in the openlo ^ ^ ^ 

y ng at the stable dynamic pressure of 75 psf is 


approximately 0.025. The weights for each operating point were scaled to obtain a unity LQR cost for a 6 
inch/sec RMS white noise disturbance input. 


Table 3.2. Cost function weights. Grayed operating points used only for fault tolerant design. 


Operating Point 


State Weight (Q | ) 

Control Weight (Q 2 ) 

Mach 0.50 

75 psf 

diag[ 1.2x1 O' 2 

I.2xl0* 2 12 12] 

610 

Mach0.S0 

132 psf 

diag[5.0xl0* 3 

5.0xl0* 3 3.5 33] 

500 

Mach 0.50 

150 psf 

diag[5.0xl0* 3 

S.QxlO" 3 2.5 2.5] 

750 

Mach 0.50 

175 psf 

diag[4.5xl0* 3 

4.5xl0* 3 1.4 1.4] 

900 

Mach 0.50 

200 psf 

diag[5.8xl0- 3 

5.8xl0- 3 0.58 0.58] 

1754 

Mach 0.50 

225 psf 

diag[9.6xl0‘ 4 

9.6X10- 4 9.6x1 0' 2 9.6xl0' 2 ] 

4800 

Mach 0.70 

125 psf 

diag[1.3xlO* 2 

1 .3xl0* 2 6.4 6.4) 

3900 

Mach 0.70 

175 psf 

diag[ 1.9x1 0‘ 3 

1.9xl0' 3 0.56 0.56] 

5600 

Mach 0.78 

75 psf 

diag[8.8xl0* 2 

8.8x10*2 44 441 

8800 

Mach 0.78 

225 psf 

diagPJxlO- 4 

3.3X10* 4 1.6x10*2 1.6 x10* 2 1 

26000 


3.4.2. Selecting the Compensator Structure and Sample Rate 

Traditionally, the design ot a multirate compensator structure begins with a successive loop closures 
structure and then incorporates cross feed between the loops as necessary. As discussed in Section 2.4.1, 
multirate successive loop closures is best applied to problems in which the closed-loop system dynamics can be 
separated into some fast dynamics and some slow dynamics. The BACT wing however does not exhibit those 
closed-loop characteristics. Closed-loop bode plots, from control input to accelerometer outputs of the BACT 
wing in feedback with a LQ Regulator, are shown in Fig. 3.8. The LQ Regulator was designed using the cost 
function weights for the mach 0.50 75 psf operating point specified in Table 3.2. Therefore the bode plot \ 
representative of the closed-loop dynamics we are trying to achieve with the flutter suppression system. Notice 
that the closed-loop dynamics have only one peak - that associated with the pitch and plunge modes - and do not 
exhibit the fast and slow dynamics traditionally associated with successive loop closures. Consequently, a 
traditional multirate successive loop closure structure is not directly applicable to this problem. 

Instead of basing our multirate compensator structure on the closed-loop dynamics of the system, we 
selected compensator structures which used different sampling schedules to reduce either the number of 
computations or the hardware required to implement the compensator. We designed four compensators: a 
single-rate (SR), a multirate successive loop closures type (MRSLC); a multirate with multiplexed inputs 
(MRMI), and a single-rate fault tolerant (SRFT). All of these compensators are second order except the fault 
tolerant design which is fourth order. 

3.4.2. /. Single-Rate (SR) 

The single-rate compensator was designed for comparison with the other compensators. A block diagram 
of this compensator is shown in Fig. 3.9. The sample/update rate for this compensator is 50 Hz. This rate is 
approximately 10 times the frequency of the dominant pitch and plunge modes. The compensator includes a 
0.02 second computational delay, which satisfies NASA’s computational delay requirement. This was achieved 
by constraining the compensator’s direct feedthrough term to be zero. 

The state space structure of the compensator is 
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Figure 3.9 Block diagram and corresponding sampling schedule for the SR compensator 
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. and CS Cmd is the command output to the zero-order-hold, a., bi and c are the free „ a ;„ 

mE^n. ( C 3 C 2 T ntS) ° ptimized lhe olher « ains wcn constrained to the values shown. The structure 5 

is a minimum realization of the second order compensator. See [Berg Mason & Yana 19911 for a 
discussion of minimum reason, Tl* samp.ing schedule for Eqn. ,3.2) ,s shown in Rg 3^ 

J.4.2.2. Multirate Successive Loop Closures ( MRSLC ) 
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Figure 3.10. Block diagram and corresponding sampling schedule for the MRSLC compensator 



Ffctio o# Fast to Slow Sample Rales 
Figure 3.11. Computational savings with the MRSLC design 

The net result of this two loop configuration is a compensator structure just like the single-rate design 
except that the digital processor needs to update one of the digital processor states only every fourth 
sample/update period. A block diagram ot this compensator along with a diagram of its sampling schedule is 
shown in Fig. 3.10. Note that this diagram only illustrates the structure of the compensator - it is not a 
thematic of how the compensator would be implemented. When actually implemented, this compensator will 
use the same number ot D/A and A/D converters as the SR compensator, but will require 37% fewer real-time 
multiplications per unit time. 

The choice of sample/update rates for the slow loop was arbitrary within the constraints of the GMCLS. 
Our goal was simply to reduce the number of multiplications required by the compensator without significantly 
degrading its performance. The 12.5 Hz sample/update rate was chosen because it is a good compromise 
between the total number of multiplications saved by utilizing this multirate structure and the ratio of the fast to 
slow sampling rates. Figure 3.1 1 shows the percent reduction in the number ot multiplication by using the 
MRSLC design over the SR design. There is a decreasing return in computauonal savings as the ratio of the 
last to slow sampling rate increases. In the limit, the compensator degenerates to a first order compensator with 
a reduction in multiplications of 50%. Based on Fig. 3. 1 1 we chose a sampling rate rauo of 4. 

The state space structure of the compensator which was used for the optimization is 

f;i<m,n + l)l [ af p ^ lfTE AcceKm.nil 
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Figure 3. 12. Block diagram and corresponding sampling schedule for the MRMI co mpensa t o r 


CS Cmdy(m,n)l ’ey 
CSCmd,(m,n)J 0 


0 ’ fzj(m,/i)l 

c j Jlz2 ("».«)/ 


(3.3b) 


CS Cmd(m.n) = CS Cmd l/(m,n) + CS Cmd ,(/n,n) 


(3.3c) 


where z i andz 2 are the digital processor states; TE Accel and LE Accel are the acceleration inputs from the 
A/D converters; and CS Cmd is the command output to the zero-order-hold. a„ b t , and c, are the free gains 
which were optimized. The other gains were constrained to the values shown. The structure in Eqn. (3.3) 
corresponds to the successive loop closures structure of Fig. 3.10. The intermediate outputs CS Cmpy and 
CS Cmp, were added to ensure that Eqn. (3.3) corresponded to Fig. 3.10. 

J.4.2.3. Multirate with Multiplexed Inputs (MRMI) 

The multirate compensator with multiplex inputs was designed to reduce the number of A/D converters 
required to implement the SR design. In this design, the compensator state and output updates occur at 50 Hz. 
The outputs of the TE and LE accelerometers are sampled at 25 Hz with a 0.02 second delay between the 
sampling of the TE accelerometer output and the LE accelerometer output. Thus, the MRMI requires only one 
A/D converter to sample both accelerometer outputs because it can be multinlexed between the two signals. In 
addition, the digital processor gains for the MRMI compensator are periodic Jly time-varying. One set of gains 
is used when the TE accelerometer output is rumpled and another set is used when the LE accelerometer output 
is sampled. Just as in the single-rate design, the direct feedthrough terms were constrained to be zero, resulting 
in a 0.02 second computational delay. This compensator requires the same number of multiplicauons per unit 

time as the SR design but it uses only one D/A converter. Figure 3.12 shows a block diagram of the MRMI 
compensator. 

The state space structure ot the MRMI compensator is 


?|(m,n + 1)1 j" 0 
? 2 (m,n + l)j [a,( n ) 


1 ' 
a 2 (n) 


?l(m,n)l 

?2(rn.n)J 



b\{n)' 

M») 


TE Accei(m,/i) 
LE AcceKm, 


">l 


(3.4a) 


CS Cmd( nt,n) * |cj (n) 


Cl 


i(n)] 


Z\(m,n) 

Z2(m,ni 


(3.4b) 


where z\ and ? 2 are the digital processor states; TE Accel and LE Accel are the acceleration inputs from the 
A/D converters; and CS Cmd is the command output to the zero-order-hold. a,<n). b,(n). and c,(n) are the free 
gains which were optimized. These gains are functions of n because they are periodically time-varying, e.g. 

o,(n) = n,(n+2) The other gams were constrained to the values shown. The sampling schedule for Eqn. (3.4) is 
shown in Fig. 3.12. 
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Figure 3.13. Block diagram and corresponding sampling schedule for the SRFT compensator 


3.4.2. 4. Single-Rate Fault Tolerant (SRFT) I 

The single-rate fault tolerant compensator was designed to highlight the multiple plant capability of our 
synthesis algorithm. This compensator is fourth order with a sample/update rate of 200 Hz and a 0.005 second 
computational delay. A block diagram of the compensator and its corresponding sampling schedule are shown 
in Fig. 3.13. The state space representations of the SRFT compensator is similar to the 2 nd order single-rate 
compensator with the exception that the digital processor is fourth order. ! 

The SRFT compensator is fault tolerant in the sense that it stabilizes ah the plant conditions even with one 
of the accelerometers disconnected. To achieve fault tolerance for all 24 plant conditions, we optimized the 
compensator for 22 simultaneous plant conditions - as opposed to just six for the preceding designs. These J 

include the six operating points used in the previous designs evaluated at three cases each: 1 ) both TE and LE 6 

sensors active: 2) only the TE sensor active; and 3) only the LE sensor active. In addition to those 18. we 
added four more operating points at mach 0.50 evaluated for the case where only the LE sensor is active. These 
operating points are grayed in Table 3.2. 

v 

3.4.3. Designing a Stabilizing Compensator 

We used the synthesis algorithm presented in Section 2.4 to optimize the gains of the four compensators 
discussed in Section 3.4.2. The algorithm requires an initial guess for the compensator's digital processor gains 
tor which the closed-loop system, the BACT wing and compensator, is stable. The difficulty in finding these ( 

gains is that the closed-loop system must be stable at all operating points used in the optimization. > 

To get a stabilizing guess for the wing at all operating points we used a boot-strapping technique. First we 
found values of the processor gains which stabilized the BACT wing at one operating point. Then we optimized 
the gains for the wing at that one operating point using large values for the plant disturbance noise and sensor 
notse intensities. The large value of noise intensities introduced uncertainty into the plant. Consequently, the 
resulting compensator was more robust than a compensator optimized for a plant with no noise. This new set of 
processor gams always stabilized the wing at the original operating point plus at least one other operating point. 

We then used the new processor gains as the initial guess to the problem with the wing at two (or morel 
operating points. The procedure was continued unu! the compensator stabilized the plant at all the operating 
points and the problem could be solved using realistic noise intensities. 

Before beginning the bootstrapping procedure we needed to find a set of processor gains which stabilized 
the closed-loop system for at least one operating point. This was straightforward for the SR and MRSLC 
compensators. We designed a first order single-rate compensator with pitch acceleration input, CS command 
output and a sampling rate of 50 Hz. Recall from Section 3.2.2 that pitch acceleration is essentially the 
difference in the TE and LE accelerations. The pole location and gain value of this compensator were found 
using root locus. The initial stabilizing guess for the SR design consisted of this first order compensator in 
parallel with an arbitrary first order compensator that had an input/output gain of zero. For the MRSLC system. 


I 
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LT a ^> f,m r 0rder " 0n,penSat ° r a$ 3,1 initial guess for *e fast loop of the successive loop closures structure, 
and an arbitrary first order compensator, with an input/output gain of zero, for the slow loop. 

An initial guess for the MRMI processor gains was more difficult to find than for the SR and MRSLC 

17~ , Due *° ,tS C ° mp,ex sampHng schedule could not design an initial guess by traditional 
' nS * ’ WC esigne<1 a compensator with ,he multiplexed structure but with very small gains Then 
we used die bootstrapping technique, beginning with the BACT wing operating at a low dynamic presT" 

and col^r l ** COmpcnsator gains were ver y small, they did not destabilize the wing 

vemes ole , f “ Z ‘"L ^ bo ‘*‘ rapping for this compensator took several iterations, 

erses one or two for the other compensators, because we began with such a poor initial guess 

To obtain an initial guess for the SRFT processor gains we began by designing two 2 »d order 

H C Stab,liZed H thC P,an * WhC " 1,16 LE «« disconnected, the other stabilized the plant 
.. t . , - Was ,sconnect ed. We then combined these two compensators into a single 4 th order design and 

wh^nlt 6 * gainS . UnU, . lhe nCW fourth order compensator stabilized the plant when both sensors were active or 
when only one or the other was active. Finally this design was used in the bootstrapping procedure discussed 
ear ter to obtain a single fourth order compensator which stabilized the wing at all operating points. 

3.4.4. Optimizing the Digital Processor Gains 

SectiTnlT U TT, Ze<1 ! he dlgUal PrOCCSS ° r ga,ns of the three compensators with the algorithm discussed in 
Section 2.4. The opumization used the following parameters: 

Plant Conditions. S,x simultaneous operating points for the second order designs; 22 

simultaneous operating points for the fourth order design. See Table 3 2 and 
Section 3.4.2.4 

The second order designs used the cost function weights listed in Table 3.2. 
The fourth order design used the weights in Table 3.2 for cases where both 
the TE and LE sensors were active, and one-tenth those values for cases 
where either sensor was inactive 

36 in 2 /sec 2 - this is the intensity of the white noise input to the Drvden filter 
and was specified by NASA 

0 rad 2 /sec 4 for initial designs. 240 rad 2 /sec 4 for final designs. This is 
discrete sensor noise for the TE and LE acceleration measurements 

Obtained using root locus and boc strapping, see Section 3.4.3 
See equations (3.2M3.4) 

See Figures 3.9, 3. 10, 3.12 and 3.13. 

In all designs the direct feed through terms were const? .lined to be zero. 

Additional gain constraints for each compensator are specified in 
Section 3.4.2. 

■ bove *■ ■" « « 

JAl Deign Itcratka Bated on Pcrfonnaac* and R^huUmss Analysis 

metiTH Symhe 7 n V he mUl,irate Compensa,ors we evaluate ‘ J ‘heir performance and robustness using the 
methods discussed ,n Section 2.5. One of the robustness measures was the maximum singular value of the 

|mnimum destabilizing multiplicative uncertainty at the compensator inputs (a structured singular value). When 
we synthesized the compensators using a sensor noise covariance intensity of zero, the size of the destabilizing 


Cost Function Weights: 


Process Noise PSD value: 

Sensor Noise PSD value: 

Initial Stabilizing Gains: 
Compensator Structure: 
Sampling Schedule: 

Gain Constraints: 
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gain was unacceptably small - less than 0.20 for the B ACT wing at some operating points. NASA had specified 
a value of 0.75. To improve the robustness at the compensator input we increased the sensor noise intensity to 
240 rad 2 /sec 4 and re-optimized the processor gains. This procedure was motivated by the Loop Transfer 
Recover technique for LQG systems described in [Doyle & Stein 1981]. The results of increasing the sensor 
noise are discussed in the following Section. 

3.5. DESIGN RESULTS 

We designed four compensators using the approach discussed in the previous sections. For review, the four 
are the: 

1 ) Single-Rate 2 nd Order (SR) 

2) Multirate 2 nd Order Successive Loop Closures (MRSLC) 

3) Multirate 2 nd Order Multiplexed Input (MRMI) 

4) Single-Rate Fault Tolerant (SRFT) 

The structure of each of these compensators was discussed in Section 3.4.2. Optimum values for the digital 
processor gains are given in Appendix A. 

We looked at five performance and robustness measures: 

1 ) Cost function value 

2) Gust pulse response 

3) Maximum RMS gain from disturbance to the control surface deflection and deflection rate 

4) Gain and phase margins at the compensator output 

5) The maximum singular value of the minimum destabilizing multiplicative uncertainty at the 
compensator input 

Results are presented for three operating points, mach 0.50 132 psf, mach 0.70 146 psf, and mach 0.78 
151 psf. Each of these operating points is 5 psf above the critical flutter dynamic pressure for the corresponding 
mach number, and so the BACT wing is nominally unstable at each of these operating points. It is important to 
note that none of these operating points were used for the compensator optimization. Therefore the 
compensators were not tuned to these particular operating points. In general, the performance and robustness of 
the compensators at these three operating points is indicative of their performance at the remaining 21 operating 
points. 


3.5.1. Cost Function Value 

One measure of the overall steady-state performance of a compensator is the value of the cost function in 
Eqn (3. 1 ) at the optimum value of the digital processor gains. (A value for the cost function is returned by our 
synthesis algorithm at the completion of the optimization.) For our 2 nd order designs, a “perfect” compensator 
would have a cost function value of 6, assuming no sensor noise. The “perfect” fault tolerant design would 
have a cost of 7.6 since it optimizes a different cost function. By “perfect” compensators we mean continuous 
LQR designs with gain scheduling, i.e., they use a different set of feedback gains at every operating point. We 
expect the costs associated with our compensators to be higher since they used discrete sampling, did not use 
gain scheduling, and had fictitious sensor noise. 

It is more realistic to compare the cost of our compensators to that of a discrete LQG design with fictitious 
sensor noise and gain scheduling. This comparison eliminates some of the differences due to sampling and 
fictitious sensor noise. The cost associated with the discrete LQG compensator is the lowest cost we can expect 
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screle LQG and tor our lour designs The cons assoc, aied with our second order compensator are alntosi 

a totT LQG J ' S,8n ™ S ' S “ *• <*— LOO ,s significantly more 

complex - it is a 16 ,fl order compensator with gain scheduling. 

3.5.2. The Gust Pulse Response 

The gust pulse response provides an indication of the transient response of the closed-loop system due to a 
.. S !! r 3 " Ce m P ut- 8 ust P u * se response was found by simulating the response of the BACT wing in 
feedback w.ththe flutter suppression system to a disturbance input pulse with an amplitude of 10 in/sec and a 
duration o, 0 004 second, Ms simuhwo. was perform* using , he M-ffl. mnta djenhd in Ailacl^nU 
Figures 3.14-3.16 show the response of the BACT wing at mach 0.70 and 146 psf to the specified 
disturbance gust pulse. Also shown is the response of the wing with a continuous LQ regulator. The cost 
unction weights for this LQ regulator design satisfy the same design criterion as was used to optimize the 
compensators gain, ,S«« Seclion 3.4.,., We proved response plols tor only one opemimg pJTg^ 
pulse responses ai other operating points are similar 10 those provided in Figs. 3. 14-3. |6. 

For comparison we also provided a gust pulse response plot lor the 2»d order compensators synthesned 

rob T CU "h' 1S ReCa " ll “‘ fic, “ io " s " 0I ” ad d' d lo the sensors in order to improve die 

46 nTd“ “ COmP '" Sa, ° r , " pu, F '^‘ 3,7 sh ”» s «* P"o" response of die BACT wing a, mach 0.70 and 
146 pst due to a gus. pulse dislurtance. The pnmtdy effect of adding sensor noise is to decrease die damping 

m0d ' S “ m °" PreVale, “ “ p,Kh res F° nse «i» in the 

The gust pulse response plots are shown below. 
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Figure 3.18. Steady-state covariance propagation at mach 0.70 146 psf with 1 in/sec RMS white noise disturbance 
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Figure 3. 19. Block diagram of discrete system for calculating RMS gain and corresponding sampling schedule 


3 . 59 . RMS Gain for Control Surface Deflection and Deflection Rate 

One of NASA's specifications was . nit on the steady-state covariance of the control surface deflection 
and deflection rate for a 1 in/sec RMS white noise disturbance. Our closed-loop system consists of a continuous 
plant and a discrete compensator. Therefore these steady-state covariances are periodically time-varying. In 
Fig. 3.18 we show the steady-state covariance propagation for the BACT wing in feedback with the three 
compensators at an operaung point of mach 0.70 and 146 psf for a unity RMS white noise disturbance. 

We calculated the values of the steady-state covariance at the sample/update times using the method 
described in Section 2.5.4. Between the sample/update times of the compensator, the covariances were 
propagated using the dynamics of the open-loop continuous BACT wing. The steady-state covariances are only 
shown for one BTP of the compensator • they repeat themselves during every BTP of the compensator. 

One meaningful interpretation of NASA's specification would be to look at the peak steady-state 
covariance value taken from this covanance plot. This value, though, is an upper limit on the closed-loop gain 
for a white noise disturbance and is not an accurate indicator of the control activity level. A better measure of 
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control activity would be the maximum RMS gam calculated using Eqn. 2.25. This is an exact measure of the 
maximum RMS gain tor any non-decaying input signal. 

In order to apply Eqn. (2.25). which is lor a discrete system, to our mixed continuous/discrete system we 
created a new discrete multirate system in which the continuous inputs and outputs of interest are sampled very 
last (see Section 2.5.3). We chose a sampling rate for the CS deflection and deflection rate of 1000 Hz. This is 
more than twenty times the control surface actuator rolloff frequency. A block diagram of this new discrete- 
time system, with the single-rate compensator of Eqn. (3.2), is shown in Fig. 3.19 along with its sampling 
schedule. This new system is now multirate even though the compensator is single-rate. The ETIS for this 
system has a sample/update rate of 1000 Hz and an N of 20. 

We used this new ETIS system to calculate the maximum RMS gam of the original system between the 
disturbance and the CS deflection and between the disturbance and the CS deflection rate. The maximum RMS 

gams for the BACT wing at three operating points are summarized in Table 3.3. See also the related work of 
[Sivashankar & Khargonekar 1991). 

3.5.4. Gain and Phase Margins at the Compensator Output 

Gam and phase margins were calculated at the compensator output using the ETIS and a multiloop Nyquist 
diagram. The ETIS of the plant and compensator were computed independently and then combined in series to 
torm an ETIS loop transfer function. Gain and phase margins were subsequently measured directly off the 
multiloop Nyquist plot of this function. These are traditional gain and phase margins, and assume that the gain 
and phase do not vary simultaneously. The details of this technique are given in Section 2.5.1, [Mason 1992], 
and [Mason & Berg 1992] (Attachment 2). 

The gam and phase margins for the BACT wing at three operating points are presented in Table 3.3. These 
values are typical of the margins at all 24 operating points, although the margins tend to be better at lower 
dynamic pressures and slightly worse at higher dynamic pressures. A representative Nyquist diagram is shown 

in Fig. 3.20. This particular Nyquist plot has two encirclements of the -1 point because the open-loop plant has 
two unstable poles. 
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Figure 3.21. Uncertainty Model 


Figure 3.22 ETIS uncertainty output feedback model 


3.5.5. Robustness at the Compensator Input 

The uncenainty at the compensator input was assumed to be a multiplicative perturbation of the form 
ihown in Fig. 3.21. where k \ and k 2 are complex gams. We transformed this system into die output feedback 
form traditionally used in robustness analysis using simple block diagram algebra. However when the 
compensator ,s muitirate we must use the ETIS of the plant, compensator and uncertainty. A block diagram of 
this closed-loop ETIS for the multirate flutter suppression system is shown in Fig. 3.22. G E is the loop tnusfer 
tuncuon consisting of the compensator and plant ETIS transfer functions connected in series. 

Now. given the system in the form shown in Fig. 3.22. we can calculate an exact value for the size of the 
smallest destabilizing perturbation (Doyle 1982] . First rewrite A £ in Fig. 3.22 as 

A £ = /] *1 + hk 2 (3 5) 

where l\ = diag{ 10 10... 10} with IN diagonal elements, and where h has a similar form. Then it can be 
shown that 


^ (A nun^-^ s up max p{[/, + /, e j6 jH E (e J<> )^ j for0< d£itand0< 9S2k\ (3.6) 

where 0(A min> represents the maximum magnitude of the smallest destabilizing * 1 or Jb; p , s the spectral 
radius; and H E {z N ) = (/ - C £ (z N ))~ l G E (z N ). “ ^ 

We are guaranteed that the system in Fig. 3.21 will remain stable as long as 




^^min ) 


(3.7) 


We are also guaranteed that when Eqn. (3.7) is violated, there exist values of k\ and k-> that destabilize the 
system in Fig. 3.21. 

Equation (3.6) 1$ straightforward to solve with a two dimensional search in 0 and ft The results are given 
tn Table 3.3. For comparison, the corresponding results for the design without the fictitious sensor noise are 
also given in Table 3.3 Notice that the addition of the fictitious no.se increases the maximum singular value of 
the smallest destabilizing uncenainty by as much as 60%. 

Even with the fictitious sensor noise, the robustness at the compensator inputs does not meet NASA's 
'pecification for a maximum singular value of 0.75. We could have improved the robustness at the 
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4. CONCLUSIONS AND RECOMMENDATIONS 


4.1. CONCLUSIONS 

The principle advantage of multirate control is that it gives the designer freedom to choose a sampling 
schedule which best utilizes the available hardware and software. In the flutter suppression system design, for 
example, we developed multirate controllers that provide performance comparable to a single-rate design, yet 
require either fewer real time multiplications per unit time to implement or require fewer A/D converters. 

The disadvantage of multirate control is that this additional flexibility substantially increases the 
complexity of design and analysis over the single-rate case. Undoubtedly, the lack of good design and analysis 
tools has discouraged many from applying multirate control even when the situation may be ideal for a multirate 
design. 

In this report we addressed the difficulty of multirate design and analysis by presenting a multirate design 
methodology. The methodology specifies a design approach and provides specific tools necessary to apply the 
approach to a practical problem. The tools are for modeling a multirate system, for synthesizing a multirate 
compensator which is robust to plant perturbations, and for analyzing the performance and robustness of a 
multirate system. The resulting methodology is powerful and straightforward to apply. 

To demonstrate the methodology we applied it to design several multirate compensators for NASA's 
BACT wing. Those compensators satisfy the specified design specifications and illustrate some of the benefits 
of multirate control. 

4.2. RECOMMENDATIONS FOR FUTURE RESEARCH 

1) Our synthesis algorithm currently requires a stabilizing initial guess for the digital processor gains. 
Obtaining a stabilizing initial guess for those gams can be difficult, especially when the multiple plant 
conditions capability of the algorithm is used, because the initial guess must stabilize all plant 
conditions simultaneously. Eliminating this requirement would substantially improve the algorithm's 
versatility. 

2) The singular value analysis of multirate systems leads directly to a structured singular value problem 
with repeated blocks. Calculating an exact solution to this problem is difficult for all but the simple 
tow parameter case. This is an area which needs further research. 
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APPENDIX A. DESIGN RESULTS 

Following are the state space matrices for the optimized flutter suppression system digital processors 
discussed in Section 3.0. 

A.l. SINGLE-RATE 2 nd ORDER 

STP=BTP=0.02 sec: A- 1 . See Section 3.4.2. 1 for a description of the sampling schedule. 
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z 2 (m,0) 
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1 — 0. 94601 LE Accel(m.O) 


A.2. MULTIRATE SUCCESSIVE LOOP CLOSURES 

STP=0.02 sec; BTP=0.08 sec; ,V=4. See Section 3.4.2.2 for a descnption of the sampling schedule. 

Update during first STP of the BTP: 
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Update during second STP of the BTP: 


: f (m, 2 ) = 0. 75673;, (m. 1 ) + [- 1 0. 37644 ]| 
CS Cmd(m.l) = 1 0 -4 [— 2. 5338 2.35354] 

Update during third STP of the BTP: 
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Update during fourth STP of the BTP: 


Zf(m+ 1,0) = 0. 75673? / (m,3) + [-1 0. 37644]| 
CS Cmd(m,3) = 10~ 4 [-2.5338 2.35354]! 
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Z s (m, 0 ) 


the BTP aSSUmed ‘ hat * S Updated dUFing ^ f,rst STP of BTP - but “ « can be updated during any STP of 


A3. MULTIRATE MULTIPLEXED INPUT 

STP=0.02 sec: BTP=0.04 sec; N=2. See Section 3 4.2.3 for a description of the sampling schedule 

Update during first STP of the BTP: Only the TE Accelerometer Is sampled. The LE Accel value is held 
from the previous STP. 
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Update during second STP of the BTP: 
from the previous STP. 


Only the LE Accelerometer is sampled.The TE Accel value is held 
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A.4. SINGLE-RATE FAULT TOLERANT 

STP=0.005 sec; BTP=0.005 sec; N =\ . See Section 3.4.2.4 for a description of the sampling schedule. 
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APPENDIX B. M-FILES USED TO DEFINE THE FLUTTEk 
SUPPRESSION SYSTEM SYNTHESIS PROBLEM 


B.l. PAPAABCD 

Format: [am. bm. cm. dm. vm]=PAPAabcd(&uune, rolloff. form) 

Description. Creates state space matrices defining the PAPA wing at operating point specified in fume such 
that 


x - anu + bm u 
v = cnu + dmu 


where v 


r plunge 
pitch 

plunge rate 
pitch rate 
TE accelerometer 
LE accelerometer 
command to actuator 
CS control surface 
CS control surface rate 
CS control surface accel 
mode 1 
mode 2 
mode 3 
mode 4 






J 


and 


u 


{ CS command 
Dryden filter input 


fname text variable containing the name of the operating point of interest, e.g. freon_m5__q75\ 
faame must have the same name as the file which contains the data 

rolloff frequency in rad/sec of first order anti-aliasing roll-off at the sensors. The filter has the 
form 


y filtered J+ro || 0 ff y unfiltered 
form indicates the desired form 

if form = 0 : am. bm . cm, dm is unchanged fre m original data 
I : am bm, cm. dm is block diagonal 

2 : am. bm. cm. dm is block diagonal with scaled states and outputs 
Outputs : am. bm.cm.ikn state space description of the plant 

vm transformation matrix used to obtain modal form 


48 


B.2. FSSCOMP 

Format: 


Inputs : 


— • — - * *» 

ctype specifies the desired compensator 

' fC,5 ’ P '* -Lc- IT. “ d ' scnp "°" ° f ,h ' ^ S «« k - R “ des,g„ 

FSScm '> re,un > s • description of the Mult, rate Successive Loop 
Closure design K 

mrm ' th des^g n SCmP retUmS 3 deSCr,ptl ° n ° f * he Mult,rate w/ Multiplexed Input 

‘srft then FSScmp returns a description of the Single-Rate Fault Tolerant 
design 


Outputs : 


W,!:;: * de “ np "°" of ,he com Pensa,o, a* by d» synthesis algorithm. 


BJ. VIROPT_SR OR MRMI 

Format. mropt_srORmrmi 

rnulTk^ 2 "''7 d “ S, " 8le - ra,e “"’’““'o' “ —• co| npensa,or w,,h 

*' SP and the MRMI dej 
none 


Description. 


Inputs : 
Outputs : 


Attachment * l0ba ' V ‘"' abi ' s — -V «*—*- — - and deftned in Section 3.3 


B.4. MROPT_MRSLC 
Format: mropt.mrslc 


Description: 
Inputs : 
Outputs : 


Defies , he input data for the ntulnnu. contpe.sa.or w„h success, ve loop closure f„„„. 
none 

V ‘'' iW ' S “” d b> - deHned in Secton 3.3 of 


B-5. VIROPTSRFT 

Format: 

mropt_srft 

Description: 

Defines the input data for the single-rate 

Inputs: 

none 

Outputs : 

Outputs to global variables used by 
Attachment 4 
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Reduced-Order Multirate Compensator Synthesis 


Gregory S. Mason and Martin C. Berg 
University of Washington, Seattle, Washington 98195 


. A * **** tyMbeeizing reduced-order maltirate compensator* it ptMtlid. Necessary conditions for which 
, !****“*“* p,r “* , * T **••«* “ tanmte timt quadrat* cost function are derived. An algoritkm 

f or hading compeataior parameter ealoct which satisfy the accessary condition. is detcribcd. This algorithm is 
thea used to design several ftp position coatroMers for a iwo-iiak robot ana. 


Introduction 

I N many cases, a multirate compensator can provide better 
performance than a single-rate compensator requiring the 
same number of real-time computations. Berg, for example, 
was able to reduce the steady-state rms response of states and 
controls to a disturbance for a simple mass-spring-mass system 
nearly 20V. by using a multirate compensator over a single-rate 
compensator. 1 Numerous other examples have been provided 
in the literature by Berg. 13 Amit. 4 - 5 and Yang. 6 Although 
multirate compensators can provide improved performance 
over single-rate compensators, they are also, in general, more 
complicated to design. 

The complexity of multirate compensators steins from the 
fact that they are by nature time varying, periodically time 
varying for most practical applications. Not only must design- 
ers choose multiple sampling/update rates for the compensa- 
tor, but they must also determine the parameter values for a 
time-varying compensator. 

One method for designing mu'tirate compensators is multi- 
rate linear quadratic Gaussian (LQG). 4 Multirate LQG is the 
multirate equivalent of single-rate LQG and is straightforward 
to solve because the equations governing the solution are sim- 
ilar to those for the single-rate case. Multiratc LQG. however, 
results in a full-order compensator which has periodically 
time-varying gains. For many applications full-order, time- 
varying compensators are not practical. 

A generalized algorithm for multirate synthesis (GAMS) 6 
was developed by Y ang to overcome many of the shortcomings 


of multirate LQG. Yang's algorithm can synthesize reduced- 
order multirate compensators with or without time-varying 
gains by using a numerical gradient-type search to find opti- 
mum compensator parameter values. His algorithm uses a 
finite time cost function in its problem formulation, unlike 
multirate and single-rate LQG which use an infinite time cost 
function. By using a finite time cost function. Yang's algo- 
rithm eliminates the numerical problem that arises when a 
destabilizing compensator is encountered during the numerical 
search. Even though Yang's algorithm uses a closed-form ex- 
pression for the gradient, the calculations necessary to per- 
form the gradient-type search are extremely cumbersome. 

In this paper, we present a new algorithm for synthesizing 
reduced-order multirate compensators with or without time- 
varying gains. The algorithm utilizes the compensator struc- 
ture of Yang’s algorithm, but the problem is formulated using 
an infinite time, instead of a finite time, cost function. This 
allows us to derive necessary conditions for which the multi- 
rate compensator minimizes the cost function. The equations 
for the necessary conditions are fairly simple and can be solved 
directly using a standard nonlinear equation solver, eliminat- 
ing many of the numerical complexities of Yang's algorithm. 


General Multi rate Compensator 

Before deriving the equations governing a reduced-order 
multirate compensator, we will first present the structure for a 
general multirate compensator. We restrict our discussion for 
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now to compensators with time-invariant gains and sampling/ 
update rates whose ratios are rational numbers. 

A general multirate compensator is shown in Fig. 1. Each 
input .v, output u, and state ; is sampled/ updated at a rate 
which, in general, represents the desired bandwidth of the 
input or output with which it is associated. The variable v is the 
value of > currently available to the digital processor from the 
zero-order hold: while u is the current output from the digital 
processor which is held with a zero-order hold to form the 
output u. When the sampling/ update rates have ratios which 
are rational numbers, the sampling/ update schedule is period- 
ically time varying. We define the greatest common divisor of 
all of the sampling / update periods as the shortest time period 
(STP) and the least common multiple of all of the sampling/ 
update periods as the basic time period IBTP ) (see Fig. 2) 


MASON AND BERG MtLTIRATE COMPENSATOR SYNTHESIS 

Fast Sampling Scheme t^-ltl itl itl |f| |t 
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*here u is a hold stale used to mode! the sampler and zero-or- 
der hold between u and u. The s. , . s t .„ , and s„ * are switching 
matrices for y . c , and u , respectively, that model the system's 
sampling/ update activity at the start of the *th STP. Also. 
s* (r has the form 
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0 
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0 0 


0 : 
0 I 
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where 


I if theyth '*• or w) is sampled/ updated 

at the start of the k th STP 

0 otherwise 
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Fig . 2 Example of a multirale sampling scheme. 

A more complete discussion of this compensator structure can 
be tound in Refs. 6 and 7. 

Equations (1) and (2) can be written more compactly 


M • l ~ V* ^ > | 

= D k \\ 


as 

(3) 

<4) 


where 


Equations (3) and (4) form a single-rate periodically time- 
varying system with a sampling rate of one STP and a period 
of on tBTP If N = BTP/STP , then A k = A k . ,, fl* = B k . N , 
C* and 0* =23** v . 

Even though A k , B k , C* , and D k are periodically time vary- 
ing, the multirate compensator gains, A % 5, C, ?.nd D , are time 
invariant. The periodic:ty of the multirate compensator is due 
to multirate sampling/ updating, not the compensator gains. In 
the remainder of this section, we will demonstrate how the 
time-invariant compensator gains. A % A, t\ and £>, can be 
separated from the periodic compensator matrices .4* , B k , C* , 
and D k . 

Define the composite compensator matrix as 


D A* C 
‘ l«* A k 


and factor P H as follows: 

p t ~S u PS lt -5,. 


where 


S,. 


t) t 

B A 


J. . 0 

0 s ; . 
0 0 
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s, t 0 /-j. t Ol 

I 0 / 0 
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0 / -5; , 0 

i.» o /-J., 
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( 10 ) 


Equation (6) is a key result. It allows us to factor the time- 
invariant compensator gains, the unknown parameters we 
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will solve for in the next section, out of the time-varying 
compensator. 

It is important to note the difference between P k and P in 
Eq. (6). P k (with a subscript) is a periodically time-varying 
matrix defined by Eq. (5). It includes all of the information 
about the compensator gains and the sampling/ update sched- 
ule. Pis a constant matrix which contains only the gains for the 
compensator. P k can be written in terms of Pand S u , S 2k . and 
S ik using Eq. (6). S u , S 2 *, and S 3 * are periodically time-vary- 
ing matrices which contain a description of the sampling/uD- 
date scheme. K 


It is important to keep in mind that P k in Eq. ^9) correspono* 
to the P k in Eq. (5). a periodically time-varying matrix which 
contains all of the information about the multirate compensa- 
tor gains and sampling/update rates. 

The closed-loop system is 

x* + \ = F ck x k ■+■ *1* (20) 

where 

F ck =F + GP k H (2i) 


Derivation of the Necessary Conditions 

In this section, we will use the results of the previous section 
to derive the necessary conditions for the reduced-order multi- 
rate compensator. The multirate problem to be solved is as 
follows. 

Given: 

the discretized plant model 


G ck * W + GP k H (22) 

The state covariance propagation for this system obeys 

X k + > * F ck X k Fj k + G ck RG c r k (23) 

where 


X k ♦ | * Fx k + <5 U k + 


(11) 


y k » Rk k + v k 


( 12 ) 


where P, (5 , If, and ft are obtained by discretizing the analog 
plant matrices at one STP ; w,. and v* arc discrete-time Gaus- 
sian white noise inputs; U is the control input from the com- 
pensator; and k is the sampled sensor output 
Find: 

the multirate control law vrith a prescribed dynamic order and 
sampling schedule, of the form of Eqs. (1) and (2) which 
minimizes a quadratic cost function of the form 



where E is the expected value operator, and the summation 
from 1 to N accounts for the fact that the closed-loop system 
is periodically time varying. A prescribed sampling schedule 
implies that the values of j* s y k and s u k are known. 

using Eqs. (3) and (4), it is easy to see that this problem is 
essentially a time-varying feedback problem— a time-invariant 
plant with a periodically time-varying compensator. One thing 
that makes this problem difficult is thai the compensator has 
an explicit form, that of Eqs. ( 1 ) and (2), in which only certain 
parameters. A % B, C, and Z>, can be adjusted to minimize J. 

To solve the multirate control problem, we cast it into out- 
put feedback form and follow a derivation similar to Muk- 
hopadhyay’s for the single rate case. 1 9 Using Eqs. (3) and (4) 
and Eqs. ( II ) and ( 12), we write the output feedback equations 



Equations (14-16) can be written more compactly as 

♦ Gu k ♦ Wft k (17) 

y* ■ nmi -► Vn k ( is) 

W. (19) 


X k mE[x k xi j. R =£jW] 

Equations (20-22) represent a periodically time-varying sys- 
tem with a period of one BTP . We can generate a single-rate 
system by repeated application of Eq. (20) over one BTP. 10 
The single-rate system can be written as 


X * + N * F \ftgX k G bk 


(24) 


where 


Ff>* * l)^rU + A/-2)£cU* /V- J) F (k (25) 

Gm * + F cik . ,,G f * I 

^ + ' j I G C ( k /v - | ) j 

(26) 


* 


■ - i 

Vk* i 

i 

= j 

Vk ♦ N - I I 


This single-rate system has exactly the same values for x as 
the periodically time-varying, closed-loop system at each BTP. 
However, the values of x at the intermediate STP are lost 
because x is incremented by N in Eq. (24) but only by I in 
Eq. (20). There are N such single-rate systems associated with 
Eq. (20). They can be written as 


4*.,. for i* 1,2, ,.V (27) 

^ A* >t stable, then the periodically time-varyirg system 
Eq. (20) is sts 'le." We can calculate the steady-state covari- 
ance for x using the following Lyapunov equations: 

X* - FhtXtFj, * G„/f,C£ , for * » 1,2, ,.V (28) 

r rt 0 O' 

|0 B 01 


1^0 0 0 /? | 

Note that i. periodic, that Is it varies within one BTP, but 
from BTPto BTP X, mX t , N . Once we have calculated X, at 
any k using Eq. (28), we can use Eq. (23) to propagate it over 
the BTP. This eliminates the need to solve Eq. (28) /V times. 

Now, using Eqs. (23) and (13), and the properties of the 
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operator, we can write the cost function for the 
stabilized system as (see Refs. 8 and 9) 

\ 

J = «?, Tr ( + MP>H + ( MP t H) T + {P k H) T Q 2 P k H\ X k 
~lP'V) T QzP t VR} (29) 

4 Cons,raints E 0- (23) to the cost J using 
Lagrange multipliers. A* , to obtain * 

\ 

J = ,?. Tr 1 + MP>H + (MP * H ) T +iPkH) T Q 1 P k H}X t 

* ( P. V) T QzP \ VR * A[. , '[F ek X k Fj k + G ck RGT k - X k . , J j 

(30) 

with *,«**.,. 

Necessary conditions for minimum J are 


mumy ddit ' 0n ' d must ** Positive definite for a mini- 

Substituting Eq (30) into Eq. (31) and replacing P k with 
i * “j; « +Sj* from Eq. (6), we obtain 

dj 

= 0 = Q, + MP„H + (MP k H) T 

* lP>H)T Qi P k H + Fi A k . tFtk - A* (32 ) 

for A: = 1,2 V with A* ■A*,*. 

<JA*., = ° = F '* X * F « + G ckRC/ k -X k . t (33) 

tor * = 1.2 V with X„ =X k . N . 

d'j s 

dP * 0 = 2 I) Si i* [ [Q: + C r A* . i C] />» [//A’*// r + K#fK rJ 

*• j Af r + G 7 A» . i F j A» // r j S[ k (34) 

Equations (32-34) are a set of coupled matrix equations 
They make up necessary conditions for P, which is comprised 
of the multirate compensator gain matrices A , B, C, and t) in 
Eqv^ ( I ) and (2), to minimize the cost function J. Values oM 
B, C .and t). found by solving Eqs. (32-34), can be substituted 

•-Jr <l,and <2) ’ along wi,h ,he defin i‘ion of the sampling 
schedule. s ; , . s „ , , and s , . to form the complete hme-vary* 
ing multiratc compensator. y 

To ensure that the compensator gains satisfying Eqs. (32- 
34) minimize J, we should also check that the Hessian of ) 
with respect to P is positive definite. Our present algorithm 
does not calculate the Hessian explicitly, but uses an approxt- 

«► 

Equations (32-34) were derived assuming time-invariant 
compensator gains. We can easily derive the correspond^ 
equations for periodically time-varying gains. Let * 

Am t>mb k (3J) 

with the restriction that p _ a- 

and £.-%*■£>, Define the composite periodically time-vary- 
ing compensator matrix 7 

* *, A, J ( 36 ) 


Then replace fi with c 

respect to P k to obtain * Eq <30> and wi,h 

II r fr 

dP k '°“ S| ‘H02 + C%. 1 G]/. 1 [ MVtW r + ^ K rj 

-(«'-C'A..,f]AT.Hrj s/i , . , 2 N (J?) 

Thus, for every new set of comoensatnr ...» 
new equation of the form of Eq 1 *^). * we ob,a ‘ n one 

Equations (32-34) are very similar to ih. . i 
tions. In fact, if we set S„, s,». and hS^ht" e<,U *' 

spond to a single-rate system, and N*l C ° fre ' 

results derived by Mukhopadhyay lor the single-rate £ 

Implementation 

To find a reduced-order multirate compensator that mini 
mizes the cost function J, we need to solve Eqs. (32-34) fw the 

dete™?ne,°h 8i “ nS ^ A n ° wchart of ,he algorithm used to 
determine the compensator gains is shown in Fig. 3. (Jsina the 

prescribed sampling schedule the algorithm first discretizes the 
analog plant model, analog cost function, and analog process 
no.se model. (See Ref. 2 for a discussion of the relevant 
cremation procedures.) Equations (32-34) are then solved for 
the compensator gains using a gradient-type search in Mat- 

Disoeuze the analog plant weighting mainces 
*nd process noise covariance 

X- — * — 1 

Build the matrices f.G. w,H. V.m&P 
using ( 14 ) -( 16 ) 

r * 

Get a stabilising compensator 


Calculate A, using (28) then propagate 
Jf] using (23) to obtain Xj. X,. . . . . ^ 

-- l 

Calculate using (A.4) then propagate 
A n using (A.3) to obtain , a 


Calculate using (34) 

= * 

Calculate the nep direction and length 


Calculate the neat guess for the competuator 


Stabilising 

compensator 1 


Reduce Lhe 
step si tt 


If tha Hassian it not poatuve 
dahnita notify tht asar 


X. I. f . and 0 found 


a! 
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Fit- 4 PI»Mr two-link robot am. 

T, 


-KH 


Compensator 


Two Link Ann 


Fig. 5 TLA plant/ compensator cot figuration. 


lab. - We chose a gradient-type search to solve Eqs. (32-34) 
because it allows us to easily add constraints on the parameters 
values— simple equality constraints were used to find the opti- 
mized compensators in the next section. The equations neces- 
sary to solve for the Lagrange multipliers are located in the 
Appendix, Eqs. (A3) and (A4). To ensure that the solution 
represents a minimum /, the algorithm checks that the Hessian 
of 7 with respect to the free parameters in P is positive definite 
at the solution point. 

Because Eqs. (32-34) are not valid when the closed-loop 
system is unstable, the algorithm I) must be provided with an 
initial stabilizing compensator, and 2) must result in a stabiliz- 
ing compensator at every iteration. From our experience, find- 
ing an initial stabilizing compensator is generally not a prob- 
syslems suitable for multirate control can be 
stabi ized using successive loop closure with minimal cross 
coupling between the control loops. A stabilizing multirate 
compensator can then be obtained by discretizing the individ- 
ual continuous control loops at the desired sampling rates 
When there are no constraints on its structure, a stabilizing 
compensator can also be obtained using the boot straoDina 

hod of B.»r » Fo, difficult multiratecontrofptotu 
lems. where a stabilizing compensator cannot be found using 
either ol the preceeding two methods, one can always use 
'ang * algorithm to find a stabilizing compensator and then 
switch to our algorithm to complete the optimisation. In our 
experience, Yang's algorithm usually converges to a stabilizing 
solution quickly— it is the optimization of the compensator 
parameters that is time consuming. 

To avoid the problem of destabilizing compensators during 
the iteration process, we included a check in the algorithm 
which systematically reduces the step size to ensure that the 
compensator is stabilizing. Because the gradient of the cost 
function with respect to the compensator parameters becomes 
very large near the stability boundary, the algorithm is always 
forced away from a destabilizing solution as long as it never 
steps over the stability boundary into an unstable region. 

Even though our algorithm was programmed as an inter- 
preted Matlab M-File we found that it still performed better 
.ian Yang s algorithm which runs as compiled Fonran. The 
primarv difference between the two algorithms is in the com- 
plexity ol the expression lor the gradient of J with respect to 


Table I SaamHag/update 
nitf for TLA 


7: 


Sample/ update rate, s 


8 

& 

Ti 

r 2 


0.225 

0.028125 

0.225 

0.028125 


the compensator parameters. Calculation of the gradient 
expression for Yang's problem involves diagonalization of 
tne closed-loop system and evaluation of several matrix cqua- 
tions with nested summations. Compare Eqs. (32-34) with 
Eqs. (1 J2-1 15) in Ref. 3 to see the difference in the complexity 
of the two gradient expressions. 

Two-LInlt Robot Ann Example 

a us * d a n \ a,hematic » 1 model of a planar two-link robot 
arm (TLA) to demonstrate the capabilities of our algorithm 
This is the same model used by Yang. 4 and so we were able to 
verify our results by direct comparison. A diagram of the TLA 
is shown in Fig. 4. 

The goal of our design was to control the tip position 6 of 
the arm via a multirate compensator. We used the following 

analog cost function and process noise covariance matrices 
rrom Kef. 6. 


J = lim£ 


r 

0.21 0 0 ol 


0 0 0 0 

< X T 



0 0 18.5 0 

s. ! 

0 0 0 0 


- 


-t- u 


r f 0.01 0 

L 0 0.69444 


*4 


(38) 


where 


••(3 


£fww r l 


0.69444 0 

0 0.01 


(39) 


We assumed perfect measurement and that plant distur- 
bances enter the system coincident with the confoi torques 
rne sampling/update rates are given in Table I 
Five difference compensators were designed: an analog 
LQ}R, a multirate lead/lead, an optimized multirate lead/lead 
an optimized multirate general second order, and an optimized 
single-rate general second order. We used a smooth step input 
to i,,i and t rt i defined as follows: 


1 




10.005 


1 -“■($]“' 


0.001 m. 


tsT e 

r f « 0.125 s 


(40) 


\ 


4 
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and the servo configuration shown c - ■ * 

performance of the different com F ' 8 5 to measurc the 
the TLA for the five compenw°or^?h° rS The c r ” pons * of 
The analog LOR 0r$ 15 shown *n Figs. 6a -6c. 

provided this 

sible using the cost function w£«hl£? f he rwponse P° s ’ 
The muitirate lead "eid wl fi«nH m *' r,CeS ofE< < (38) 
closures. We designed the controU^* succes$ive *o°P 
so that the eigenvalues of the ctoS d,SCWe d ° m * ,n 

.te. »« ob<«,nM l Q r ;rz,::\T 

This compensator consists of two simole i«h i„„ screte ,,me - 
j » r, opmiini .be Tj’Z 

"on, . ,o r, opm„„ „ ,b« “°~ 

The final three compensators were svnthMt 7 Mf 
new algorithm and the cost-weighting matrices used^od^gn 
the anaiog LQR compensator. The optimized multirate lead/ 
lead was found by optimizing the pole/zero locations and 

closures. thC Cad/Ie * d compcnsa,or foun <* by successive looj 

The optimized multirate general second-order compensator 
uses the same sampling/update scheme as the lead/lead com- 
pensators but has the compensator structure of Eq. ( 41 ) where 
a,, . b „ , c „ , and d„ are the parameters which were optimized 
This compensator has the maximum number of independent 
free parameters possible for a second-order system.’ 


(41) 


The optimized single-rate general second-order compensa- 
or is a single-rate equivalent of the multirate general second- 

o^^° mPe, li $at0 I ,l haS the S * me struc,Uf e as the multirate 
general second-order compensator, Eq. (41), but uses a single 

sampling rate. This sampling rate was chosen such that the 
number of computations required to implement either the mul- 

are *th °s«ne * rat * COmpensaIors durin * real-time operation 

Our results are the same as those obtained using Yang's 
algorithm. They demonstrate how multirate compensators can 
rKKrl?* i be ,,er Performance than single-rate commutators ™ 
trading lower bandwidth control of the slow modes for higher 

^ dW,dt Ji COnt ? 1 ° f ,he f “‘ modes ,n « h » example, we were 
able to reduce the tip response overshoot 404b and the pe ak 
control torque 25 V# by using a multirate controller over a 
single-rate controller. 


r «n 

L o 

0 

°22 j 

J . 6 = 

r i 
1*2. 

r] 

fc„ 

c, i 

t fi * 

U. 

rf.t] 

L 



u. 

dn\ 


In this paper, we have presented a new algorithm for synthe- 
sizing reduced-order multirate compensators. It can be used to 
design compensators of arbitrary structure and dynamic order 
with independent sampling/update rates for the compensator 

°“ t P uts ', and $,ate * Thi * algorithm provides the ver- 
tility of Yang s algorithm without the numerical complexi- 
ties associated with the finite time cost function. 

do t0 discount Yang's algorithm alto- 

gether because, while our algorithm requires an initial stabiliz- 
ing compensator, Yang’s does not. For those problems where 
finding an initial stabilizing compensator is difficult, we can 

US'.,* “** Y *Cf ‘ • l,omhm 10 find • stabilizing compensator 
and then quickly optimize the compensator parameter values 
with our algorithm. 


Wh l ch ,ubil ’“* ,he naultirate system, we can 

,h< ,l<ad '{' ,,a, « va,ue * of A ‘ Wh we A, is defined by 
Eq. (32) rewritten here as Eq. (A I). 


°* Cl MP.H * i\fP,H) r * iP t H) r Q } P,H 
+ - A* 


(Al) 


*r- 
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for k * 1,2, ,.V with A* = A* * * . 

First simplify Eq. (Al) by defining 

n . r c?i ^ 1 f r / 1 

Q ' W ftj- y * “ L IA2 > 

where / is an identity matrix. Then Eq. (Al) can be written as 

A* = Jl Q>J k + Fj k A *♦,/>* (A3) 

for k = 1,2, with A* * A***. 

Equation (A3) represents a periodically time- varying Lya- 
punov equation. We can create an equivalent single-rate sys- 
tem by repeated application of Eq. (Ai). 

= ^dk Q Jdk + Fjk ^k^dk (A4) 

for k = 1.2 .V with A*»A*„*. 

Fdk = Fctk + s- i)F cii + N ~2)F c{ i + j si ~ i) F ck (A5) 


Jd* = i 


j Ak + N-\ } Fcik + N-2)F c{k + N- J) • ■ ■ F ck "j 

•^♦♦*-2>^cU + N-3> F ck 


I ft 0 

0 Q } 


0 o, 


Equation (A4) is a time-invariant Lyapunov equation which 
can be solved for A* . Once any A* has been found, the propa- 
gation Eq. (A3) can be used to find the remaining A*. 
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Introduction 

^TUMEROUS approaches have been taken to analyze the stabil- 
M*-- robustness of multirate systems. Most notable is the 
work of Thompson.' Thompson and Dailey.* Meyer and Bumis. 3 
and Khareonekar et at 4 References i and 2 analyzed the gain and 
phase margins or a multirate system using an approach based on 
Kranc vector switch decomposiuon. Meyer and Bumis stuped the 
stability of mulurate systems and then frequency domain responses 
by applying the concept of block processing to multirate systems. 
hJiargonekar et al. analyzed the robustness of penodic compensa- 
tors. a super set of multirate compensators, using isomorphisms. 

Although these approaches seem quite diverse, fundamentally 
they are very similar. These three approaches, along with most mui- 
tirate analysis techniques, use a transformation, such as Kranc vec- 
tors or block processing, to convert the multinte system into an 
equivalent single-rate system which can be analyzed with estab- 
tshed single-rate techniques. The single-rate results are then used to 
characterize the stability and robustness of the original multirate 
system. 

* n this paper we compile the important multirate analysis results 
from Refs. 1-6 and present them in a unifying state-space formula- 
tion. In addition, we provide some new results that clarify the rela- 
tionship between a mulurate system and its single-rate equivalent. 
Finally, we apply all of these results to a pracucal example: a mulu- 
rate flutter suppression svstem designed for a model wing. 

Summary of Multirate Analysis Tools 

In this section we summarize some important and useful multi- 
rate robustness analysis results. These results are applicable to mul- 
tirate systems that are linear, causal, finite-dimensional, and whose 
sampleAipdate/delay activities are periodic and synchronized to a 
common clock. 

Remit I 

A mulurate system can be modeled as an equivalent single-rate 
system (ESRS). Modeling a multirate system as an ESRS is funda- 
mental to mulurate robustness analysis. The ESRS allows one to 
manipulate arid anal y ze a mulurate system as if it were single rate, 
using the ESRS. single-rate and multirate systems can be combined 
msenes or in feedback loops just as in classical control. 4 It has also 

° Cen * $,n ?*** r >te/mulnrate system will be stable when- 

ever its ESRS is stable. 


anS^SaSS 9I ’ 28 ' 2 * ** A,AA 
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The ESRS of a multinte system can be obtained by modeling 
the multinte system as a periodically time-varying system*-*and 
transforming the periodically ume-varying system into an M RS 3 
The state-space representation of the ESRS is 

.ttm-t-l.O) = AgAm.Oi+B^im.O) (U) 

>’f (m. 0) = CfXim. 0) + D E u t tm. 0) (J5) 

where 


y(m.O) 

y £ (m.O) = . - v(m ’ ,) u £ (m. 0) 


u(m, 0) 
u[m, 1) 


( 2 ) 


Ly (m. A/- 1 )J 


|_tt(m, N - 1) 


form « 0, 1. 2 and n = 0. 1, 2 N- \ 

In these equations, u{m, n) and y(m, n) represent the values of 
the input and output, respectively, of the original multirate system 
at the (mN+n )th sampling instant. The integer N is the ratio of the 
least common multiple of all of the multirate system’s sample/ 
update/delay periods to their greatest common divisor. The sub- 
script E denotes vectors and matrices strictly with the 

ESRS. 

A key feature of an ESRS is that us tnput/output vectors art 
composite vectors containing the tnput/output values of the multi- 
rate system at N separate sampling times. Consequently, an ESKS 
is always multiple input, multiple output (MIMO) even if the orig- 
inal multirate system is single input, single output (SISO) Another 
key feature is that an ESRS always has a nonzero direct feed 
through term D E . This is because D t contains information about 
how past inputs affect the current output. For systems with no 
dynamics, the direct feed through term D E is block diagonal. For 
example, the ESRS of a constant uncertainty matrix A is 


A f * block diag[ A. A A| ( 3 > 

with N blocks. 

Refer to Refs. 3. 6. 8. and 9 for the details on modeling a multi- 
rate system as an ESRS. 


Result 1 

A discrete signal Mm. /n is related io its ESRS signal wAm. 0) 
as follows: * ^ 


Mm, 1 > •• • W^n-N + 1 )Jwf<m, 0) 

where is a switching function defined as 
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that 


IV' vU) =1 it k = -mN form = 0. 1. . . . 

= 0 otherwise (4b> 

Result 2 provides a convenient mathematical connection 
between the input/output vectors of a multirate system and its 
ESRS. The torm ot VV v given in Eq. (4a* is useful for analytic 
results whereas that given in Eq. (4b) is convenient for numencal 
simulation. 

Result 2 allows one to use single-rate analysis techniques to 
obtain results for multirate systems. For example, to compute the 
pme domain response of a multirate system to a generic input we 
can find the time domain response of the ESRS using any applica- 
ble single-rate technique, and then transform that response back to 
the original multirate system using result 2. 

Result 3 

The two-norm and ims value of a discrete signal Mm, ft) and its 
ESRS signal w^m, 0) are equal or, equivalently, 

II w(m, n >ll : = 0)\\ 2 and rm$[Mm, ft)] *rms( 0)] 

Result 3 follows directly from Eq. (2) and the definition of the 
two-norm and rms value. 10 

Result 4 

The maximum rms gain of a multirate system is given by the H- 
infinity norm ot its ESRS transfer function G^z) or, equivalently. 


sup 

ma (a) mO 


rms [y (m, ft) ] 
rms (u(m, ft) ) 


i|Gf(*>!L 


It is well known that the maximum rms gain of a SISO single- 
rate sy?*em is equivalent to the maximum gain on that system’s 
Bode plot. Although a transfer function for a multirate system does 
not exist in the traditional sense, we can see from result 4 that the 
//-infinity norm of a SISO multirate system plays the same role as 
the maximum Bode plot gain of a single-rate system. 

Result 4 follows from result 3 and Refs! 10-12. It is also 
derived, in pan. in Ref. 4. 

Result 5 

The singular values ot a single-rate transfer function G<:) and its 
ESRS transfer function Gf<:) are related as follows: 

= | a(G(^)] r a(G(e'‘* 


where <j| | denotes a column vector of singular values. 

Result 5 relates the singular values of a single-rate system to the 
Mnguiar values of its ESRS and provides some insight into the sig- 
nificance of the singular values associated with the ESRS of a mul- 
ti™ 6 system. From result 5 we can see that one effect of trans- 
forming G(.*) into Gfi:) is that the singular values of G<:) at high 
frequencies ;ire biased into G,(:i at lower freauencies. Conse- 
quently. the u\m. m in result 4 that results in the maximum rms 
gain docs not necessarily contain the frequency u> associated with 
; G r is the sampling pc nod of the ESRS). The input sig- 

nal of maximum rms gain must be constructed using the right sin- 
gular vectors of G^;) and result 2. We will demonstrate this proce- 
dure in ihe following section. (See the Appendix for a denvation of 
result 5.) 

Result i 

The stability, gam marg in. and phase margin of a SISO multirate 

■ 'item *.ait be determined dirccth from j N .quist plot of its ESRS 
Recall from result l that the ESRS ot a SISO multirate system \s 
MIMO. therefore, me muttiloop iNyquist staoility criterion must 


be used in result 6. 13 This, however, is one case where gain and 
phase margins taken from a muiuloop Nyquist plot can be inter- 
preted in the traditional sense because gain and phase variations at 
the multirate system's input/output apply simultaneously to ail the 
inputs/outputs of the ESRS. 

Result 6 follows from Eq. (3) and is derived in Ref. 1. 


Result 7 

The robustness of a multirate system can be determined by 
applying structured and unstructured singular value analysis to that 
system’s ESRS. Given the ESRS transfer function G e of the nomi- 
nal system and the uncertainty transfer function A* we can apply 
established singular value analysis techniques to find the sine of 
the smallest uncertainty <*(A f ) that destah iii™^ the closed-loop 
system in Fig. 1. This result however, is only a conservative esti- 
mate of the size of the actual smallest dgnta*Ji»ywi| uncertainly A. 
The input/output vectors of an ESRS are composite vectors, con- 
taining the input/output values of the multirate system at N tuple 
times. Thus <r ( Ag) found using unstructured singular value analy- 
sis accounts for not only the fictitious perturbation normally asso- 
ciated with unstructured singular values, but also far ttme-vmy mg 
and noncausal perturbations. A valid pert ur bation far a given 
d(A*) might, for example, include block diagonal cta me ptt m A* 
that are unequal. This corresponds to a time-varying 
because the gain between u(m,/i) and y(m, a) varies with a. 
Another valid perturbation could include no nzer o upper b lock 
diagonal elements in Af. This corresponds to a noncausal perturbs 
don because a future input u(m, n+l)can affect the c u rren t output 
y(m, a). 

For the ESRS uncertainty A c to represent the ***** » uncertainty 
A, its structure must obey Eq. (1). Finding fr(A*) subject to Eq. 
( 1 ) requires the solution of a structured singular value problem. 
Unfortunately, even simply structured dynamic unesnainties in a 
singte-ratc/multirate system transform to uncertainties with com- 
plex structures in the ESRS. The complex structure makes it diffi- 
cult to obtain a good estimate of the size of the destabiliz- 

ing structured perturbation. However, when the single-rate/ 
multirate uncertainty is a constant, as is the case for many prob- 
lems, the ESRS uncertainty is also a constant with a repe a le d block 
diagonal form (see Eq. (3)J. A good estimate of the solution of 
such a structured singular value problem with repeated blocks can 
be found using one of the methods in Refs. 14-17. 

Result 7 follows directly from the fact that a single-rate/multi- 
rate system is stable if and only if its ESRS is stable. 7 

Application 

In this section we apply the results of the previous section to a 
real world example: a multirate flutter suppression system for a 
model wing. This application points out some of the practicalities 
of muiurate robustness analysis. 

The model wing used in this example is being developed under 
the Benchmark Active Controls Project at the NASA Lmgky 
Research Center. It consists of a ngid airfoil mounted on a pitch 
and plunge apparatus (PAPA). The PAPA mourn provides the two 
degrees of freedom needed to model classical wing flutter. The 
wing has one control surface located on the trailing edge (TE) of 
the airfoil. Two accelerometers measure pitch and plunge accelera- 
tions. We used a 1 5th -order mathematical model of the wing for 
the control system design. This model incorporates a second-order 
Dryden gust niter, a third-order actuator model, and a lOth-order 
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airfoil model (two rigid body modes and six unsteady aero states). 
A block diagram of this model is shown in Fig. 2. 

A second-order multirate flutter suppression system was 
designed for the model wing. The multirate compensator, shown in 
Fig. 3, consists of two hnt-order control loops. The slow loop is 
sampled/updated at 50 Hz. The fast loop is sampledAipdated at 200 
Hz, resulting in N = 4. We optimized all of the free parameters in 
this compensator using the multirate compensator synthesis algo- 
rithm described in Ref. 9. The free parameters were the pole and 
zero locations for the fast and slow loops, the gain values for the 
fast and slow loops, and gam values for the cross feed between the 
fast and slow loops. 

To analyze the robustness of the closed-loop system we exam- 
ined the gam and phase margins at the plant input and output and 
the rvns gains from disturbance input to the control surface deflec- 
tion and deflection rate. Gain and phase margins provide a measure 
of the uncertainty allowed in the plant model. ' r he nns gains pro- 
vide a measure of the allowable disturbance level before the con- 
trol surface actuator limits are exceeded. 

Tra dition a l gain and phase margins at the plant input were cal- 
culated using result 6. The location of the co r res p o nding complex 
loop gain k M is shown in Fig. 4. The multiloop Nyquist diagram for 
the open-loop ESRS is shown in Fig. 5. Gain and phase margins 
calculated from the Nyquist diagram are given in Table 1. 

Generalized gam and phase margins at the plant eutput were cal - 
culated using the ESRS and the structured singular value. For this 
analysis, the closed-loop system was cast into the form of the 
block diagram shown in Fig. 1 . The uncertainty block A and the 
corresponding A, are 


k , 01 
_° 


jV. <M 

i° vj 


(3) 


wIiwb /, it an N x N identity matrix. The complex gain uncertain* 
tiea *| and 4, rt pree ent additive plant output uncertainties. They 
an shown in Fig. 4. G t (in Fig. I ) ia the nominal closed- loop ,y». 
'em .-omonsed of the model wing and the muitirate flutter suppres- 
sion system. 


The form of the uncertainty A f in Eq. (5) leads directly to a 
structured singular value problem with repeated block 

uncertainties. Typically, an exact solution for such a problem is 
difficult to find. Fortunately, this problem has only two free param- 
eters. k, and 4., and an exact solution can be found. We r«ir.ii.t^f 
an exact value of the size of the smallest da«t»hii mn e uncertaint y 
using the following equation' 4 : 


1 


<J(A„ 


- = sup max p {(/(©) C,(e /4 ) } 

,) 0«*«« 0 < S < 2, 


(6) 


with 


t/< 0) 


!/ 0 i 
Lo Aj 


where a (A*.) is the maximum singular value of the umIIm 
d est abilizing A; and p< ) is the spectral radius. The resulting value 
of #(A^) is given in Table 1. 

Generalized gain and phaae margins co rre sp on di ng to *( A _j 
were cal cu la t ed usiug the method described in Ref. 18. Figure 6 
shows the “region of guaranteed stability” for simultaneous (inde- 
pendent) gain and phase changes in 1 +4, and 1+k,. 
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Table I Ro bmunw remlti for multir»t* nuttw suppremon ,y WMB 

{ idin margin at plant input (with *, =* : =oi, dB [~T2 ' + 10) 

Phase margin at plant input t with *. =k : =0K dee 
with 4r,= | 

Maximum rms gam. v to 6^. deg s/in. 

Maximum rms gam, v to 6^. deg/in. 


923 


r65 

0.74 

0.125 

4.67 


Themaxunom mu gain values for disturbance input to control 
Mirtace deflection and deflection rate were calculated 

' R .!f Ul1 4 ’ however * “ d ««ly applicable only to discrete 
ems and not to our mixed continuous/discrete system. There- 

deflwn^^S^li 6 PUW * i3mbmct “P« •** control surface 

outpu,s M500Hz - assuming a zero- 
order liold at the disturbance input. This sampling rate is 30 times 

rms gains are then the peak values of the maximum singular value 
plots ot this new muitirue system (see Fig. 7). The peak values, 
which occur at m^- 27 rad/s. are given in Table 1. Since we 
assumed a stair step disturbance input, the values in Table 1 are 
lower, yet very dose, bounds on the maximum rms gains. 19 

ontrary to what occurs in linear time inva riant single-rate sys- 
tems. we cannot assume that the signal producing the maximum 
gam n) has the simple form sinOv—.r/n). in«^ 

we must construct n) using result 2. For our example 

^™j WM slgTUj ° f m * ximum "ns gain n> from thedis- 

|«rt»nce^ to the control surface deflection ^ was found as fol- 

G £ (e^w r )=transfer function from v to ^ 

evaluated at <u mn (7g) 


and 


a y*'i 


«,^i 



fa this paper we have summarized some important 
robustness analysis results. The foundation for these resuks is the 
equivalent single-rate system. This system allows one to analyze 
the stability and robustness of a muitiraie system using well- 
known single-rate techniques. 1 

There are. however, drawbacks to using the equivalent single- 
rale system. First, the composite structure of its inputs end otnm 
leads to mamx transfer hincuons of high order. This is a mta 
especially v hen the ratio of the longest to the shortest snmpiiag/ 
u pdat e/delay period of the multinue system is huge. Second, even 
when the simplest constant parameter uncertainty m ode l is used 
for robustness analysis, the composite input and output str uct ure 
gives rise to a structured singular value problem with repealed 
diagonal blocks. This is a numerically difficult problem to solve 

fa short multixare system stabibty and robustness analysis baaed 

on the equivalent single-rate system is stnughtforwred the 

analysis can be performed using established single-rare (ecteiques, 
but the results must be interpreted carefully in accordnce with its 
input/output structure. 


Lemma: 


Appendix 


-nght singular vector associated with o (G c (e'“muT)| (7b) 
^ ,cn - «> is given by 


<r(G e (<-- v »>l = {(r{G<e*)] T 

o[G(e , '*** mm )] T ■ • cr(G(e /l ** ,JW, ' v * ,, ' v ")]r|r 

Proof: Let G(z) and Gf(r v ) be the transfer functions of a sin 
gle-rate system and its ESRS. respectively, such that 


- ( W w (n) Mf y (a - I ) . . . IV v (n - jV+ |)j 

otjSinfoi^Tm + a, ) j 
x i®s « B (*nmTm + * t ) I 

a v sin.o). I)ll 7>»*0 v ) J 


y(z) « GUMz) and y £ (/') * G^u^z*) (Al) 


. , / ■ •tMfcuwi vi f dociuk me rarapunc 

period for the ESRS is N times that of the single-rale system 
cuued with G(z). 

It is shown in Ref. 3 that 


**)•[/ :- J / :■ *-'t] yi{ r') (A2) 

where / is an identity matrix ot appropriate dimensions. 

Now define 


where H-'y is the switching function described in result 2. and T 
the sampling penod of the ESRS. 

From Eqs. (4i and 18) it is straightforward to see that, in genera 

max,mum «•'" for * mulnrate system is con 
pnsed ot the sum of sinusoids of several distinct frequencies. I 
this example, though. n) is comprised almost purely of 

Mngle sinusoid of frequency of 27 rad/s. This it because higher fn 
guenev signals which might interact with ihe multirate compensi 
• >r and increase the rms ram are atienuaied bv the Drvden nlte 
io puni o\njfTiiCN 
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and 


’ 1 
«(♦ :>J 


Gui* block diag|G(<b’;i. G(<t»:i 

G(4"'.)| 

(A4) 
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where 
such that 




v(r)=G(r)5(;» 


<A5) 


usmg E™ A™ Wn “ ***** “* 5( ‘*’ m ,enns 01 v * ,r ' ) 


where 


>(:) T ( z ) y £ ( r s ) 2(:)=7'(r)u f (r v ) 


(A6) 
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.. 

^ - i N - 1 1 
' * - / 



i -i 

i - (V- |> 


TC )= / 

(6 :) / 

(0 Z) / 

(A7) 

J 

( 0 V * .. 

V — l ~ < N — 1 } 

• (0 z) / 1 



Then from Eqs. ( A 1 > and (A5) through (A6> 

G f <z / '')=r- 1 <5(;:)r 


(A8) 


Noncmg that T < e *) W = NI . we can define a un.tarv matrix 
/ such that 

T (*>*)= vRT ( ti *) and 7Ye'V=<l/vW)f(e'»)* (A9) 

Now. using Eqs. (A8) and (A9) and the propcntes of umtarv 
matrices, it is straightforward to show that 

ofGjfe^)] =»o|<5 («^*)] (Al 0) 

“ block l«e E<)- (A4)]. and the singular values 

° f lu l dia(?onal matm *** *e union of the singular values of 
each block, we can rewrite Eq. < AIO) as follows: 

o(G f (e'' v *!]= |a(G(e>*)] y cr[G(e rt ** l2 ’ ,v 'l)| r 

<x|Gfe ,, *'*' ,4 * <w, )) r • . • o(G(e' l *' , ' ,Jl "' y-l,,,y ii)jr|r 
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abstract 

A new methodology for multirate control system desig.i i^ described. It 
accommodates a general multiple-input multiple-output control law structure that allows 
the sampling rates for the plant sensor output signals, the update rates for the processor 
states, and the update rates for the plant control input signals to be independently 
specified. It includes a capability to design for multiple plant conditions so as to achieve 
robustness to plant parameter variations. Its analysis components include a method for 
determining conventional gain and phase margins, a method for determining a bound on 
the smallest destabilizing uncertainty, and a method for determining the maximum RMS 
gain of a multirate system. The methodology is demonstrated by application to the 
design of a multirate flutter suppression system for a model wing. 

INTRODUCTION 

Multirate control systems occur frequently in engineering practice. They have 
received comparatively little attention in the technical literature. There are several 
reasons for this. One is a lack of recognition by the research community of the practical 
motivations for multirate controllers. Compared to single-rate controllers, they offer the 
relatively obvious real-time computing efficiencies in multi-loop, multi-function, multi - 
time-scale systems. But with real-time computing hardware costs as low as they are, such 
efficiencies usually do not justify the additional complexity of a multirate design. 
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In practice, multirate controllers are often necessitated by hardware constraints. For 
example, when a sensor provides a signal that is updated only at a fixed interval, except 
when the update period happens to be a suitable sampling period for a single-rate 
controller, a multirate controller must be used. 

The logistical burden that a multriate system presents is a second reason for the lack 
of attention given to multirate control systems. This burden is a consequence of the fact 
that a multirate system is time-varying from one sampling instant to the next. 
Fortunately, a well designed software package can spare the designer from most of the 
burden of the logistical difficulties, thereby allowing him or her to concentrate on the 
more fundamental design issues. 

This paper describes a multirate control system design methodology for which we 
have developed such a software package. The methodology was originally proposed in 
Reference 1. It employs the control law synthesis algorithm described in Reference 2, 
and the modeling and analysis tools described in References 3 and 4. The description is 
via an application to the design of a multirate flutter suppression system for a model 
wing. 

The remainder of the paper is divided into four sections. The first describes the 
model wing, its open-loop characteristics, and the flutter suppression system design 
goals. The second describes the design methodology and its application to the flutter 
suppression system design. The third presents the results of the flutter suppression 
system design. Conclusions are given the final section. 

PROBLEM DESCRIPTION 
The BACT Wing 

The Benchmark Active Controls Technology (BACT) Wing is being developed at the 
NASA Langley Research Center to study the modeling, prediction, and control of 
aerodynamic flutter. It consists of a rigid airfoil mounted on a flexible base. The base, 
called the Pitch and Plunge Apparatus (PAPA), provides the two degrees of freedom 
necessary to model classical wing flutter. The airfoil has one control surface (CS) located 
on the trailing edge. Two accelerometers, one near the leading edge (LE) and one near 
the trailing edge (TE) measure the airfoil's motion. References 5-6 describe the BACT 
Wing in detail. 

The flutter suppression system was designed using 16* order linear state models of 
the BACT Wing developed by NASA Langley's Structural Dynamics Division. Each 



model consists of 4 rigid body states corresponding to the pitch and plunge modes 
6 unsteady aerodynamic states, a second order actuator model, a second order Dryden 
filter, and two first order anti-aliasing filters. Figure 1 shows a block diagram of this 
structure. NASA provided us with 24 such models, each describing the dynamics of the 
wing in Freon at a different operating point. The operating points include dynamic 
pressures above and below the critical flutter pressure at three different mach numbers. 
See Table 1 for a summary of the operating points. 

Table 1. Openuing points for BACT Wing. All operating points assume Freon medium 

Dynamic Pressure (psf) 

(Unsta ble operating points are in gray) 

MachO50 lT* loo 122 132 150 175 200 225* 

Mach 0.70 75 100 125* 136 146 175* 200 225 

-Mach 0.78 7 5* 100 125 1 41 151 175 200 225* 

* Operating points used for compensatoi^synthesi^™^™^™** 111 " 1 ^™**’™* 8 *^™ 

Open-Loop Characteristics of the BACT Wing 

Two modes - pitch and plunge - dominate the open-loop dynamics of the BACT 
Wing model. For example, the poles and zeros of the CS command to the LE and TE 
accelerometer output transfer functions at mach 0.5 and 75 psf are shown in Figures 2(a) 
and 2(b). As dynamic pressure increases, one pair of these poles moves toward the right 
half plane and crosses the imaginary axis at the stability boundary. Figure 3 shows this 

pole movement. The corresponding movements of the open-loop poles not shown in 
Figure 3 are relatively small. 

The dominant pitch and plunge modes are observable at all operating points with 
either the TE or the LE accelerometer outputs, and are controllable at all operating points 
using the CS command input. The zeros of the CS command to TE accelerometer and the 
CS command to LE accelerometer transfer functions are shown in Figure 2 for the mach 
0.5 and 75 psf operating point. As dynamic pressure increases, the non-minimum phase 
zeros associated with the TE accelerometer migrate into the left half plane. The 
minimum phase zeros, associated with the LE accelerometer and located near the 
dominant poles, migrate into the right half plane. See Figure 3. 

, At ,OW dynamic P ressures transfer functions from the CS command input to the 
T£ and LE accelerometer outputs are non-minimum phase. Non-minimum phase 
systems are more difficult to control than minimum phase systems?. An alternative 
output is one which measures the difference between the two accelerometer signals. This 
new output is essentially pitch acceleration. The CS command to pitch acceleration 
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transfer function is minimum phase for all operating points, and the BACT Wing is 
relatively easy to control using this new output. We chose not to use this output directly 
in our designs, however, because it does not adequately account for the inevitable 
uncertainty in the TE and LE acceleration measurements. 

Design Objectives 

Our primary objective was to design a multirate flutter suppression system that will 
stabilize the BACT Wing when it is flown (at some future date) in the wind tunnel at 
speeds between mach 0.5 and 0.78 and dynamic pressures between 75 psf and 225 psf. In 
addition, the following constraints, most of which are functions of the hardware that will 
be used to implement the control law, were specified by NASA. 

Control Activity Constraint : For unity RMS white noise input disturbance ( 1 in/sec 
RMS), the steady state covariance of the CS deflection must not exceed 0.0625 

deg 2 (0.25 deg RMS), and the CS deflection rate must not exceed 65 deg 2 /sec 2 
(8.0 deg/sec RMS). 

Sampling Rate Restrictions: The minimum sampling period is 0.005 sec. For 
multirate sampling, all sampling periods must be multiples of 0.005 sec. 

Computational Delay : All compensators must account for a minimum 0.005 sec 
computational delay. 

Robustness Constraints: The gain and phase margins at the compensator output, 
which is a scalar signal, must be at least ±6db and ±45°. The maximum singular 
value of the smallest destabilizing multiplicative uncertainty at the compensator 
input must be 0.75, which corresponds to simultaneous gain and phase margins at 
the two sensor inputs to the compensator of ±6db and ±45° 8 . 

Finally, our multirate control law was to provide the same performance and stability 
robustness as a comparable single-rate controller yet require less hardware to implement. 

The Design Methodology 

We designed two flutter suppression systems for the BACT Wing: a single-rate 
system, for use as a baseline for comparison, and a multirate system. Each was designed 
using the methodology in Reference 1. This methodology defines a general approach and 
provides the specific tools needed to solve a design problem. The methodology has three 
parts: modeling the multirate system, optimizing the digital processor gains, and 
analyzing the performance and robustness of the closed-loop system. In the following 

paragraphs we describe the methodology and its application to the design of the two 
flutter suppression systems. 
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Modeling the Flutter Suppression System 

A block diagram of a generic flutter suppression system is shown in Figure 1. In this 
system, each sampler at the plant output can operate at an independent rate, the digital 
processor can update each processor state at an independent rate, and each zerc-order- 
hold at the compensator output can operate at an independent rate. This compensator 
model has the form of the Generalized Multirate Control Law Structure (GMCLS) 
discussed in References 1, 3 and 4. The GMCLS provides a framework for modeling 
multirate compensators and eliminates much of the bookkeeping involved with multiple 
sample/update rates and/or time delays. Using the GMCLS, it is straightforward to 
represent a multirate system as a periodically time-varying single-rate system. Later we 
will see that the resulting periodically time-varying system can be further transformed 
into a time invariant system by “lifting” 9 or “block processing” 10 . 

There are two components to the GMCLS: the “sampling schedule” and the “digital 
processor gains ’. The “sampling schedule” indicates the sequence of sample and update 
activities for all samplers, processor states, and zero-order- holds. In the GMCLS, all 
sample and update activities must occur at integer multiples of a specified time period T\ 
and the sampling schedule must repeat itself every NT, where N is an integer. Often T 
and N are functions of the hardware used to implement the control law. The second 
component of the GMCLS is the “digital processor gains” A z , B z , C z , andD,. These gains 
determine the dynamics of the digital processor and are typically free design parameters. 

We modeled both the single-rate and the multirate flutter suppression systems using 
the framework of the GMCLS. Both compensators have the form of the generic 
compensator shown in Figure 1, with two inputs, TE and LE accelerations, one output, 
CS command, and second order digital processor dynamics. 

The single-rate compensator has a sample/update rate of 50 Hz, which is 
approximately 10 times the frequency of the dominant pitch and plunge modes. The state 
space structure of the compensator’s digital processor is 


w»+m ro i]f •*,(«)! 

\x 2 (n+l )J ~|a, 

CScmd W-h 


+ 


0 b\ fTE Accel(n)1 

1 bi |LE Accel(n)} 


(la) 

(lb) 


where x\ and x 2 are the digital processor states; TE Accel and LE Accel are the 
acceleration inputs from the A/D converters; and CS cmd is the command output to the 
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zero-order-hold. a t . b iy and c, are the free gains to be optimized. The structure of (1) 
represents a minimal realization of a second order compensator^. 

The GMCLS sampling schedule corresponding to the single-rate compensator is 
shown in Figure 4. In this figure, circles on each time line indicate when a particular 
sample or update activity occurs. The GMCLS “digital processor gains” correspond to 
the matrix elements in the digital processor’s state space description given in (1). 

The multirate compensator was designed to provide the same performance and 
stability robustness as the single-rate compensator using a reduced number of analog-to- 
digital converters. In the multirate compensator, the digital processor states and CS 
output are updated at 50 Hz, while the accelerometers are sampled at 25 Hz. In addition, 
there is a 0.02 second delay between the sampling of the TE accelerometer output and the 
LE accelerometer output. Consequently, the multirate compensator requires only one 
A/D converter to sample both accelerometer outputs. 

To maximize the benefits of its multiplexed sampling schedule, the multirate 
compensator uses periodically time-varying digital processor gains. One set of gains is 
used when sampling the TE accelerometei output and another set is used when sampling 

the LE accelerometer output. The state space structure of the multirate compensator’s 
digital processor is thus 


fj^Oi+1)] 

JO 1 ' 

J* l(w) l + 

o />,( n y 

fTE Accel(n)) 

(2a) 

|x 2 («+l)J 

[aj(n) a 2 (n) 

ItaWJ 

1 l> 2 (n) 

[LE Accel(rt)} 

• CScmd(u)=[c,(n) c 2 (n)]j 

1*2 («>J 



(2b) 


where x\ and x 2 are the digital processor states; TE Accel and LE Accel are the 
acceleration inputs from the A/D converters; and CScmd is the command output to the 
zero-order-hold. a,(n), />,(*), and c,(n) are the free gains to be optimized. These gains are 
functions of n and are periodically time varying, e.g„ a,(n) = a,(n+2). 

The sampling schedule for the multirate compensator is shown in Figure 5. Notice 
the multiplexed sampling scheme of the two A/D converters. The GMCLS “digital 
processor gains” correspond to the matrix elements in the digital processor’s state space 
description given in (2), just as in the single-rate case. 

Optimizing the Digital Processor Gains 

To determine the values of the digital processor gains for the two compensators, we 
used the low order multirate compensator synthesis algorithm described in Reference. 2. 
along with the multiple plant conditions idea of Ly * ‘-'2. The synthesis algorithm uses 
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numerical optimization to determine values of the digital processor gains that minimize a 
quadratic cost function. The multiple plant conditions idea employs a cost function 
which is the sum of the costs associated with a single compensator in feedback with a 
nominal plant and perturbed variations of that plant. The digital processor gains that 
minimize this new cost function are robust in that they stabilize the nominal plant and the 
specified variations of that plant. 

The multiple plant conditions cost function has the form 



where £ is the expected value operator; the integer N p is the number of simultaneous 
plant perturbations under consideration; the vectors x, and n, represent the plant states and 
control inputs, respectively, of the ith plant condition; and Q x , M„ and /?, are the 
weighting matrices associated with the ith plant condition, and are free parameters 
selected by the designer. 

Using the synthesis algorithm in Reference 2 and the cost function in Eqn. (3), we 
found values of the digital processor gains for the single-rate and multriate flutter 
suppression systems that stabilized the BACT Wing at all 24 operating points. Instead of 
optimizing over all 24 operating points, however, we selected six representative ones. 
The six are indicated in Table 1. 

For each of the six operating points, we selected a unique set of weights, M,-, and 
ft;, for the cost function in (3). To select Q x , M x , and /?, we used that fact that the gains 
which minimize (3) corresponds to the LQ regulator feedback gains when the 
compensator is a continuous time design, the plant outputs are the values of the plant 
states, and the compensator is strictly a feedback gain. Accordingly, we first designed a 
continuous LQ regulator, one for each of the six operating points, that satisfied NASA’s 
performance criterion. The weights were chosen so that the closed-loop damping of the 
pitch and plunge modes was greater than 0.07, and the RMS control surface activity 
constraints specified by NASA were satisfied. (For comparison, the damping in the 
open-loop BACT Wing at a stable dynamic pressure of 75 psf is approximately 0.025.) 
Next, we uniformly scaled Q, M„ and /?, to obtain a unity LQ cost for each operating 
point for a 6 inch/sec RMS white noise disturbance input . Finally, we used the values of 
Qi, A/„ and R, from each LQ regulator design for the corresponding values in (3). 
Optimum values of the digital processor gains were found by minimizing the cost 
function (3) using the synthesis algorithm in Reference 2. 
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After optimization, we evaluated each compensator’s performance and robustness 
using the methods discussed in the next section. One of the robustness measures we 
considered was the maximum singular value of the minimum destabilizing multiplicative 
uncertainty at the compensator inputs (a structured singular value). The size of the 
structured singular value for our initial designs was unacceptably small. It was less than 
0.20 at some operating points, whereas NASA had specified a value of 0.75. To remedy 
this we added fictitious sensor noise and reoptimized the processor gains (our initial 
designs assumed no sensor noise). Our addition of fictitious sensor noise was motivated 
by the Loop Transfer Recover technique for LQG systems design 13 . 

Analyzing the Flutter Suppression Systems 

Recall that we modeled the flutter suppression systems using the GMCLS. It is easy 
to represent a compensator having the GMCLS as a periodically time-varying single-rate 
system. Furthermore, it is straightforward to transform a single-rate periodically time- 
varying system into a single-rate time-invariant system using “lifting” or “block 
processing . We refer to the resulting single-rate time-invariant system as the Equivalent 

Single fete System (ESRS). See References 1 or 4 for a discussion of the properties of 
the ESRS. 

An ESRS has the following form 

x(k+N)=A E x(k)+B E u E {k) (4a) 

yE(k)=CE-x(k)+D£U E (k ) ( 45 ) 

where 


/£(*) = < 

><*) ' 
M*+l) 1 

► and u E (k)=< 

f u(k) 

1 «(k+ 1) 


y{k+N-l) 


n(k+A-l)j 


and where **) are the states of the discrete system, u(k) are the discrete inputs, and \{k) 

are the sampled outputs. The subscript E denotes vectors and matrices strictly associated 
with the ESRS. 

A key feature of an ESRS is that its input/output vectors are composite vectors 
containing the input/output values of the original system at N sample times. 
Consequently, an ESRS is always MIMO even if the original system is SISO. Another 
feature is that an ESRS always has a nonzero direct feed through term. When the original 
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system has no dynamics, the direct feed through term of its ESRS. D E , is block diagonal. 
For example, the ESRS of a constant matrix A is 

Af = Block Diag[ A, A A] with N blocks (5) 

The ESRS allows one to manipulate and analyze single-rate and multirate systems as 
if they were both single-rate. ESRS state space or corresponding transfer function 
descriptions can be used to calculate input -output relations for systems in series or in 
feedback loops just as in classical control 14 . For example, to calculate the ESRS of a 
time-invariant plant in feedback with a multirate compensator, we would calculate the 
ESRS of the plant and compensator individually, using the same value of T and N for 
both, and then combine them using block diagram algebra. Furthermore, we could 
determine the stability of the original closed-loop system by calculating the eigenvalues 
of the new closed-loop ESRS system *5. 

We used the ESRS to evaluate the performance and robustness of our multirate and 
single-rate flutter suppression systems. First, we formed their closed-loop ESRS’s, and 
then we applied analysis techniques for linear time-invariant systems to the resulting 
single-rate systems. In the following section, we discuss the results of these analyses. 

DESIGN RESULTS 

By way of review, our two flutter suppression systems are: 

1) Single-Rate Second-Order with TE and LE acceleration inputs and CS command 
output 

2) Mulurate Second-Order with multiplexed TE and LE acceleration inputs and CS 
command output 

We compared the performance and robustness of these two compensators in the 
following areas 

1) Gust pulse response 

2) Maximum RMS gain from disturbance to the control surface deflection and 
deflection rate 

3) Gain and phase margins at the compensator output 

4' The maximum singular vmue >f the minimum destabilizing multiplicative 
uncertainty at the compensator input 
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The results are presented for three operating points, mach0.50 132 psf, mach0.70 
146 psf, and maeh 0.78 151 psf. Each of these is 5 psf above the critical flutter dynamic 
pressure for the corresponding mach number, so the BACT Wing is nominally unstable at 
each of these operating points. It is important to note that none of these operating points 
were used for the compensator optimization and so the compensators were not tuned to 
these particular operating points. Although we will discuss performance and robustness 
at only three operating points, the following results are indicative of the compensator’s 
performance and robustness at all 24 operating points. 

Gust Pulse Response 

The gust pulse response provides an indication of the transient response of the closed- 
loop system to a disturbance input. We computed the gust pulse response by simulating 
the response of the BACT Wing in feedback with the flutter suppression system to a 
disturbance input pulse with an amplitude of 10 in/sec and a duration of 0.004 seconds. 

Figures 6-7 show the response of the BACT Wing at mach 0.70 and 146 psf to the 
specified disturbance gust pulse. Also shown is the response of the wing with a 
continuous LQ regulator. The cost function weights for this LQ regulator design satisfy 
the same design criterion as were used in the multirate and single-rate designs. The gust 
pulse responses at the other operating points are similar to those shown in Figures 6-7. 

For comparison, gust pulse response plots for the multirate compensator synthesized 
without fictitious sensor noise are also shown in Figures 6-7. Recall that we added 
fictitious sensor noise to the multirate design in order to improve the robustness at the 
compensator input. The primary effect of adding sensor noise is to decrease the damping 
of the pitch and plunge modes. As can be seen, the reduction in damping is more 
prevalent in the pitch response than in the plunge response. 

Max RMS Gains 

One of NASA s specifications was a limit on the steady state covariance of the 
control surface deflection and deflection rate for a 1 in/sec RMS white noise disturbance. 
Our closed-loop system consists of a continuous plant and a discrete compensator. 
Therefore, these steady state covariances are periodically time varying. In Figure 8, we 
show the steady state covariance propagation for the BACT Wing in feedback with the 
two compensators at an operating point of mach 0.70 and 146 psf for a unity RMS white 
noise disturbance. 

One meaningful interpretation of NASA’s specification would be to look at the peak 
steady state covariance value taken from the covariance plot. This value, though, is an 
upper limit on the closed-loop gain for a white noise disturbance and is not a true 
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indicator of the control activity level. A better measure of control activity would be the 
maximum RMS gain. 

The maximum RMS gain of a multirate system is given by 

sup su _ RMSfrgflfc)) N 

RMSOo*oRMS(u(*)) RMS(Uf)? 60 RMS(w £ (fc)) £ Z ^°° ^ 


where IIG£(z^)ll 00 is the H-infinity norm of the transfer function of the ESRS system 

between the input and output of interest, u and y respectively. See References 1 and 4 for 
details. 

To apply the discrete equation (6) to our mixed continuous/discrete system, we 
created a new discrete multirate system in which the continuous inputs and outputs of 
interest are sampled very fast. We chose a sampling rate for the CS deflection and 
deflection rate of 500 Hz. This is more than ten times the control surface actuator roll-off 
frequency. Figure 9 shows a block diagram and sampling schedule of this new discrete 
time system in feedback with the single-rate compensator. This closed-loop system is 
multirate even though the compensator is single-rate. The ESRS for the system has a 
sample/update rate of 500 Hz and an N of 10. 

We used this new ESRS system to estimate the maximum RMS gain of the original 
single-rate system between the disturbance and the CS deflection and deflection rate. A 
similar method was used to calculate the maximum RMS gains for the disturbance input 
to the CS deflection and deflection rate for the multirate flutter c uppression system. The 

results, for the BACT Wing at the three representative operating points, are summarized 
in Table 2. 


Gain and Phase Margins 

Gain and phase margins were calculated at the compensator output using the ESRS 
and a multiloop Nyquist diagram. The ESRS of the plant and compensator were 
computed independently and then combined in series to form the ESRS loop transfer 
function. Gam and phase margins were subsequently measured directly off the multiloop 
Nyquist plot of this function. These are traditional gain and phase margins, and assume 

that the gain and phase cannot vary simultaneously. The details of this technique are 
given in References 4 and 1 6. 

The gain and phase margins for the BACT Wing at three operating points are 
presented in Table 1. These values are typical of the margins at all 24 operating points, 
although the margins tend to be better at lower dynamic pressures and slightly worse at 
higher dynamic pressures. A representative Nyquist diagram is shown in Figure 10. This 
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Table 2. Performance and Robustness Summary 



Mach 0.50, 1 32 psf 

Mach 0.70, 146 psf 

Mach 0.78, 151 psf 


Single-rate 

Multirate 

Single-rate 

Multirate 

Single-rate 

Multirate 

Max RMS Gain 
Distr. to CS Deflect. 

0.22 

0.25 

0.19 

0.19 

0.11 

0.11 

(deg sec/in) 

Max RMS Gain 
Distr. to CS Def-rate 

6.5 

6.9 

2.4 

2.6 

1.5 

1.5 

(deg/in) 







Gain Margin 

12 db 

12 db 

10 db 

10 db 

9 db 

9db 

Phase Margin 

41° 

38° 

45° 

40° 

43° 

40° 


0.41 

0.45 

0.38 

0.44 

0.35 

0.45 

c 

o 

£ 







J*1 01 

0.25 

0.35 

0.26 

0.32 

0.25 

0.31 

c 

o 

£ 







w/o sensor noise 








particular Nyquist plot has two encirclements of the - 1 point because the open-loop plant 
has two unstable poles. 


Robustness at the Compensator Input 

The uncertainty at the compensator input was assumed to be a multiplicative 
perturbation of the form shown in Figure 11, where k\ and complex gains. We 

transformed this system into the output feedback form traditionally used in robustness 
analysis using simple block diagram algebra. However, when the compensator is 
multirate we must use the ESRS of the plant, compensator and uncertainty. A block 
diagram of this closed-loop ESRS for the multirate flutter suppression system is shown in 
Figure 12. Ge is the loop transfer function comprised of the compensator and plant 
ESRS transfer functions connected in scries. 

Now, given the system in the form shown in Figure 12, we can calculate an exact 
value for the size of the smallest destabilizing perturbation 17 . First rewrite Ae in 

Figure 1 2 as 


Ae= l\k\ + hki 


(7) 
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where /i = diag{ 1 0 1 0 ... 1 0} with IN diagonal elements, and where I 2 has a similar 
form. Then it can be shown that 17 


) — 


sup max p |[/j +/2* 70 ]w£(f?^)j 


v-1 


( 8 ) 


where oiAmm) represents the maximum magnitude of the smallest destabilizing or k?, 
0 < <(> < it: 0 < 0 < 2n: p is the spectral radius; and H^z N ) = G£z N \l- 
We are guaranteed that the system in Figure 1 1 will remain stable if 

fo *°J < * (A ™ ) < 9 > 


We are also guaranteed that when (9) is violated, there exist values of k\ and ki that 
destabilize the system in Figure 1 1 . 

Equation (8) is straightforward to solve with a two dimensional search in 0 and 0. 
The results are given in Table 2. For comparison, the corresponding results for the design 
without the fictitious sensor noise are also given in Table 2. Notice that the addition of 
the fictitious noise increases the maximum singular value of the smallest destabilizing 
uncertainty by as much as 60%. 

Even with the fictitious sensor noise, the robustness at the compensator inputs does 
not meet NASA’s specification of a maximum singular value of 0.75. We could have 
improved the robustness at the compensator output further by increasing the fictitious 
sensor noise level, but we chose not to do so because doing so simultaneously reduces the 
gain and phase margins at the compensator output. 

Conclusions 

A new methodology for multirate control system design has been developed. It 
accommodates a general multiple-input multiple -output control law structure that allows 
the sampling rates for the plant sensor output signals, the update rates for the processor 
states, and the update rates tor the plant control input signals to be independently 
specified. It includes a capability to synthesize a single control law for multiple plant 
conditions so as to achieve robustness to plant parameter variations. Its analysis 
components include a method for determining conventional gain and phase margins, a 
method for determining a bound on the smallest destabilizing uncenainty, and a method 
for determining the maximum RMS gain of a multirate system. As is demonstrated in 


Mason & Berg Multirate Flutter Suppression System Design for a Model Wing 14 

this paper by application to the design of a multirate flutter suppression system for a 

model wing, this new methodology is a practical and effective tool for multirate control 
system design. 
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Figure I . Block diagram of BACT Win$ in feedback with a generic flutter suppression system 
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Figure 4. Sampling schedule for single-rate compensator 
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Figure 6. Closed-loop plunge response to gust pulse disturbance at mach 0.70. 146 psf 
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Figure 7. Closed- toop pitch response to glut pulse disturbance 
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Figure 8. Propagation of steady srate covariance of CS deflection at mach 0.70, 146 psf 
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Figure 12. ESRS for the system in Figure 1 1 
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Section one 

Introduction 


1.0 OVERVIEW 

The design and analysis of even a simple multirate compensator can be a complex task. In this document we 
escribe some Matlab M -Files which aid in multirate design. They include M-Files for modeling multirate systems 
for computing optimum values of a multirate compensator’s gains, and for generating time domain simulations 
TTte remainder of the document is divided into three sections. In Section Two we review the basics of multirate 
modeling and our optimization algorithm. We also present the key concepts and notation which are utilized in 
Section -nvree. Section Three describes each M-File. detailing its inputs and outputs. IHp M-Files in this section 

are divtded mto three categories: modeling, simulation, and synthesis. Finally, in Section Four we conclude with a 
multirate design example. 


U SOFTWARE REQUIREMENTS 

In addition to the standard Matlab toolbox routines, this software uses M-Files from both the Control_Toolbox 
and the Opumizauon.Toolbox. These two toolboxes must be present for the multirate software to operate. 


section two 

Background 


2.0 OVERVIEW 

In the following paragraphs we present the notation used in the remaining sections and review some key 
concepts which will be helpful to those using this software. A detailed explanation of the ' ... ry behind our M-Files 
can be found in References 1-3. 

2.1 A MULTIRATE FEEDBACK SYSTEM 

A multirate sampled-data system consists of a continuous plant in feedback with a multirate compensator A 
bl ck diagram of such a system is shown in Figure 2.1. The vector y s in this figure represents the continuous plant 
sensor output. Each element of y s can be sampled at an independent rate. The vector y represents the sampled 
value of y s available to the digital processor. (In our double index notation, the discrete signal p(m,n) results from 
sampling the continuous signal p(t) at the times t=(mN+ n)T ; where the integer N is the period of repetition; T is 

the sampling time; m =0, 1 , ... ; and n=0, 1 AM.) The digital processor obtains the current « aluc of y and 

combines it with the current processor state vector, f, using the state space structure shown in Figure 2.1. Each 
element of the processor state vector, ?, can be updated at an independent rate. The continuous output from the 
compensator, represented by the vector u, is formed by holding the output from the digital processor.*, with a zero- 
order-hold. Each element of the vector u can be held at an independent rate to form u. The v^tors v and w 
represent the discrete sensor and continuous process noise respectively. 

Conceptually, one can divide the multirate compensator into two parts, the “sampling schedule” and the digital 
processor gains. This is the approach used in the synthesis and analysis software. The “sampling schedule” is a 
description of when each compensator input is sampled and when each output and processor state is updated, while 

the digital processor gains determine the dynamics of the digital processor. In the following paragraphs we discuss 
each in detail. 



Figure 2. 1 Closed-loop Multirate System 
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2.1.1 The Sampling Schedule 

In general, the sampling and updating of v 5 , f , and u in Figure 2. 1 can occur at any time. We. however, require 
that these sample and update activities occur only at integer multiples of some fixed time, called the shortest time 
period (STP). The actual value of the STP is arbitrary, but is often a function of the hardware and software used to 
implement the control law. We also require that the sampling and updating activities of the sensors, states and 
outputs repeat themselves after some fixed period of time. The period of repetition of the sampling schedule is 
called the baste time period ( BTP). Finally, we define 

BTP 

the integer N = — and the value 7= STP (2 .1) 


In our double index notation, the first index <m) in p(m,n) indicates the integer number of BTPs which have elapsed 
when the sample/update occurred and the second index («) indicates the integer number of STP’s which have 
elapsed within the current BTP when the sample/update occurred. 

We can represent the sampling schedule for the multirate compensator graphically, as shown in Figure 2.2. 
e figure shows a time line for each sampler, processor state, and zero-order-hold. The time line is divided into 
one sTP increments. On the left side of the time line is a description of the signal or state, being sampled or 
updated. On the nght side is a description of the particular activity represented by the time line. e.g. state update, 
sampler, or zero-order-hold. Circles on each time line indicate when a sample or update activity associated with that 

particular signal or state takes place. Usually the sampling schedule is shown for only one BTP since the sampling 
schedule repeats itself every BTP. e 

In most applications, the sampling/updating activities for a given sensor, output or state will be periodic within 

Litf 17 ’ ** IS Sh ° Wn ,n Figure 2 2 However - **“ sampling/updating activities do not have to be periodic within the 
BTP. The only requirement is that the sampling/updating activities have some period of repetition (the BTP) and 
that they occur at integer multiples of the STP. Once the STP and BTP have been selected, the designer can 
arbitrarily specify sampling/updating activities at any mu'tiple of the STP with in one BTP. An example of a 
mulUrate sampling schedule in which the sampling/updating activities are not periodic within the BTP is shown in 

Figure 2.3. A sampling policy like this might be used to multiplex multiple inputs through a single analog to digital 
converter. © & 


2.1.2 Digital Processor Gains 

The processor gams are the values of the matrices /»-. B z . C-. and D z in Figure 2. 1 . Like the sampling schedule, 
they can be periodically time-varying with a period of repetition of one BTP. Generally, these matrices are free 
esign parameters which can be adjusted by the designer tc improve the performance of the multirate compensator 
The synthesis software discussed later in this document can be used to calculate optimum values for these gains. 

12 THE EQUIVALENT TIME-INVARIANT SYSTEM (ETIS) 

A multirate compensator with the structure discussed in Section 2.1 can be modeled as a periodically time- 

varying single-rate compensator b’ appending appropriate hold states to the digital processor model. This new 
compensator has the form 

dm,n+l) = /(|/iWwi)+ R(n)uim.n) (22a) 

y s {m,n)=C(n)x(m,n)+ D(n)u(m/u (2 .2b) 

for n = 0. 1 ... JV- 1 , and m = 0, 1 ... . 

This compensator has a sampling period of one STP and a period of repetition of one BTP (or NT), v f(mjt) 
represents the values of sit) sampled every STP; u(t) is formed by holding u(rn.n) with a zero-order-hold. 
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Figure 2.2. Example multirate sampling schedule with periodic sampling/updating activity 
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Figure 2.3. Example multirate sampling schedule with aperiodic sampling activity 


The periodically time-varying single-rate compensator can be further transformed into a single-rate Equivalent 
Time-Invariant System (ETIS) with the form shown below. 

xf m+ 1 ,0; = A e x[ m,0) + B ,0) 

>'£(m,0) = Cgt(Bi,0) + D gu/j{m,0) 

where 


and y£(m, 0) = 

Vj(/n,0) 
.Vi(m, 1) 

i/£(m,0) = 

u(m, 0) " 
u(m.l) 


_y,(m,/V-l)_ 


-u(mJV-l)- 


(2.3a) 

(2.3b) 


We use the subscript £ to denote vectors and matrices strictly associated with the ETIS. Notice that the ETIS 
input/output vectors are composite vectors containing the input/output values of the multirate (or periodically-time 
varying) system at N sampling times. Consequently, an ETIS is always MIMO even if the original system is SISO. 
If the multirate system has p inputs, q outputs and a sampling period of one STP then the ETIS is a single-rate linear 
time-invariant system with Np inputs, Nq outputs and a sampling period of one BTP. 

The ETIS is fundamental to the analysis of multirate systems. It allows one to evaluate the performance and 
stability of complex systems comprised of multirate, periodically time-varying and/or single-rate components using 
only modified linear time-invariant single-rate techniques. For example, to evaluate the stability of the system in 
Figure 2.1, we would first transform the multirate compensator into its ETIS with a given value for N. Then we 
would discretize the plant at the STP of the compensator using a zcro-ordcr-hold and transform the resulting single- 
rate system into an ETIS using the BTP of the compensator. Next, the plar* and compensator ETIS’s could be 
combined in feedback just as if they were traditional single-rate systems. Finally, we could determine the stability 
of the original multirate sampled-data system from the eigenvalues of its closed-loop ETIS. 
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23 OPTIMIZING THE DIGITAL PROCESSOR GAINS 

The synthesis software discussed later calculates optimum values of the digital processor gains, A z , B z , C\, and 
D : , by minimizing a quadratic cost function which reflects the peiiormance of the closed-loop system in Figure 2.1. 
The cost function has the form 


7- 



u T (t) 


<21 


M' 

Qi 


■*<»n 

,M(0 Jj 


(2.5) 


where 7 is the cost associated with the closed-loop system shown in Figure 2.1. The vector y c is the continuous 
criterion output and u is the continuous control input Q\ , Qi and M are the cost function weighting matrices which 
are selected by the designer. 

This cost function has the same form as that minimized by a continuous time LQR design. Thus the cost 
associated with the optimized multirate compensator and that of an LQR design can be compared directly. The 
designer can also use this fact to help select appropriate values for Q \ , Qi and M . 

To improve the robustness of the compensator, the optimization algorithm can optimize the digital processor 
gains for several different plant perturbations simultaneously. The resulting compensator will stabilize the plant at 
each perturbation and provide overall optimum performance. This is accomplished by minimizing the new cost 
function ot Eqn. (2.6) which is the sum of the costs associated with each plant perturbation. 



Here J, is the cost associated with the i th plant perturbation and there are Np plant perturbations. 
Optimum values of A z , B z , C z , and D z , occur when 


dJ 

"5a z 



= 0. 



and 


dJ 

3Z5T 


= o 


(2.7) 


Our algorithm uses a gradient type numerical search to determine values of the digital processor gains such that the 
conditions m Eqn. (2.7) are satisfied. Because the synthesis software uses an iterative process to determine optimum 
values for the digital processor gains, the user must provide the software with an initial guess for A Z ,B Z , C z , and D,. 
The compensator corresponding to these values must stabilize every plant perturbation considered in (2.6). 

Obtaining a suitable stabilizing initial guess can sometimes be a difficult problem. We refer the interested reader to 
Reference 1. 


Section three 


M-File Documentation 


3.0 STANDARD VARIABLE DEFINITIONS 

Many of the M-Files discussed in this document require similar input variables or provide similar outputs. To 
simplify the documentation of these M-Files a set of standard variables are used throughout this document. They are 
defined in Table 3.1 with Matlab variables and functions bolded. 


Variable 


Table 3.1 Standard Variable Definitions 


Description 


Plant matrices in the form pit = [F, G; //, J] where the state-space representation of the 
plant is 



i(f) = Fx(t) + Gu(t) 

(3.1a) 


y{t) = Hx(t)+ JuU) 

(3.1b) 

or 

4w+D= Fx(m,n)+ Gu(m,n ) 

(3.1c) 


yim,n) = Hx(mji)+ Ju(m,n) 

(3. Id) 


depending on the value of stp defined later. 


nplt Number of states in (3.1), or equivalently the number of rows in F 


emp 


Multirate compensator gain matrices. Given the state-space representation of the digital 
processor 

z(m,n + 1) = A z (n)z{m,n) + B z (n)y(m,n) (3.2a) 

u(m,n) = C z (n)z(m,n) + D z (n)y(nun) (3.2b) 


where m=Q,l,... and n=0,l,...^V-l, so that the gain matrices, A z , B z , C z and D z , are 
periodically time-varying, emp has the form 


cmp = 


T/UO) fl,(0)' 

'A-(l) fl-(l)' 

r 

Jc z (0) Dj(0)j 

C z ( 1) D,(l) 

’ * 


A,(N-l) B Z (N- 1) 
C(/V- 1) D.(Af-l) 


(3.3) 


If A z , B z , C v and D : are constant, then 


emp = 


A z ( 0) B z ( 0) 
C z { 0) D : (0)J 


(3.4) 


and it is assume that A z (0)=A z ( 1)= ... 1), £,(0)=£ : (1)= ... =fl,(/V- 1), etc. 

The software automatically deduces from the size of emp, the value of nemp and the 
number of compensator inputs whether the digital processor gains are periodically time- 
varying or not. 


Continued on following page ... 
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Table 3. 1 Standard Variable Definitions < continued from previous page) 

Variable Description 

ncmp Number of digital processor states in Eqn. (3.2), or equivalently the number of rows in A : . 

stp If stp > 0 then stp is STP, the shortest sample/update period of the compensator. See 

Section 2.1.1, and pit describes a continuous plant. 

If stp = 0 then the plant matrices pit describe a discrete system whose sampling rate is 
one STP of the compensator, and pit describes a discrete plant. 

stppbtp Number of STP’s per BTP. stppbtp = N. See Section 2.1.1. 


su 


sv 


sz 


Sampling schedule for tne compens ‘or output, su has two forms 

Case I: su is a matrix, su must have as many columns as there are compensator outputs, 
and must have N rows (A = BTP/STP). Each element of su has a value of 1 or 0. A 1 
in the »* column and j^ 1 row of su indicates that the i* compensator ou'out is updated 
at the STP within the current BTP. A 0 indicates no update and the previous value 
is held. For example, given a two output system with the following sampling schedule 

Output #1 — | — ? ZOH 0 

Output #2 | £ J 20H Y = update activity 

STP — J U. N=3 

— J| — BTP 

the corresponding value of su is 


su = 


1 

0 

1 


0 

1 

1 


(3.5) 


Case II: su is a row vector. The i^ 1 column of su specifies the sampling period of the 
compensator output hold in number of STP’s. For example, su=[ I 4] specifies that the 
first output is updated every STP and the second output is updated every fourth STP. 
This form assumes that the updating is synchronized on the first STP of the BTP. 


Sampling schedule for the compensator input, sy has the same form as su except that it 
specifies when the plant output (the compensator input) is sampled. 

Sampling schedule for the compensator states, sz has the same form as su except that it 
specifies when the compensator's digital processor states are updated. 



9 


3.1 MODELING M-FILES 


3.1.1 mr2etis 

Format: [ae.be ,ce ,de )=mr2etis( cmp, ncmp.su, sy, sz, stppbtp; 

Description: Convens a multirate system into an ETIS with the form: 




v(m, 0) 

x(m + 1,0) = aex(m, 0) + be^ 

y(m,n + 1) 



y(m,yV - 1) 

u(m, 0) 1 


v(m,0) 

u(/w,m + l) 

► = cex{ m, 0) + de< 

y(m % n+ I) 

u{m,N - 1) 


y(m.N - 1) 


Inputs. cmp.ncmp.su. sy . sz, stppbtp The compensator description. See Table 3.1. 

Outputs: ae.be.ce.de The ETIS matrices in Eqn. (3.6) 


3.1.2 mr2ptv 


Format: 

Description: 


[P. nu. ny, nz]=mr2ptv(cmp, ncmp.su. sy , sz. stppbtp) 

Converts a multirate system into a periodically time-varying system of the form 

jrf/w.n+l)= A(n)x(m.n ) + B(n)\\mji) 
u(m,n)= C(n)x(m,n)+ Din)y(ni.n) 


Inputs : cmp, ncmp , su, sy , sz. stppbtp The compensator description. See Table 3.1. 

Outputs: P The periodically time-varying system matrices in Eqn. (3.7), where 


To(0) c(0)l 

'£>(1) <7(1)' 


~D(N- 1) C( N - 1)1" 

_[fl(0) /t(0)j 

.fl(l) /t(l) 


B(N- 1) /t(N-l)J 


nu Number of inputs u in Eqn. (3.7) 

ny Number of outputs y in Eqn. ( 3.7) 

nz Number of states jc in Eqn. (3.7) 


(3.6a) 


(3.6b) 


(3.7a) 

(3.7a) 


(3.8) 


10 


3.1.3 ptv2etis 

Format: [ae, be, ce, det=ptv2etis(P, nu, ny , nz, stppbtp ) 

Description: Converts a periodically time-varying system into its ETIS of the form of Eqn. (3.6) 

Inpu , : P, nu. ny , nz The periodically time- varying system matrices. See Section 3.1.2 

stppbtp See Table 3.1 

Outputs . ae, be , ce, de The ETIS matrices. See Eqn. (3.6) 


3.1.4 sr2etis 


Format : 
Description : 


Inputs : 
Outputs : 


[ae, be, ce, de^=sr2etis(a, b, c, d, stppbtp ) 

Converts the single-rate system in Eqn. (3.9) into an ETIS with an N of stppbtp, and the form of 
Eqn. (3.6) 


Mm,n +1)= a(n)x(nun) + b{n)y{m,n) 
u(m,n) = dn\x(m,n)+ d(n)y(m,n) 

a. b, c, d The single-rate systems matrices. See Eqn. (3.9) 
stppbtp The desired N of the new ETIS. See Table 3. 1 

ae, be, ce, de The ETIS matrices. See Eqn. (3.6) 


(3.9a) 

(3.9b) 


3.1.5 stack 


Format: 

Description: 

Inputs : 


Outputs : 


(yh= stack (w. stppbtp) 

Converts a traditional discrete-time vector w into an ETIS vector y 

w A matrix whose columns contain the values of the vector y(nun) for Nk successive samples times, 
where A: is an integer. The matrix w has the form 


w = [y(0,0), v(0, 1 ), ... ,V(0,AM) y(l,0), y( 1,1) yt\.N-\) ><*-l.AM)) 

>' ' s written using the double index notation of Section 2.1. 
stppbtp See Table 3.1. 

y An ETIS vector of the values of v£ (the ETIS of y) with the form 

y(0.0) 

v(O.l) 

y(0.A/-l)J 



y( 1,0) ‘ 


y(k - 1.0) ' 



y( 1.1) 


y(k - 1, 1) 



xi./V-D 


! v(Ar — 1, N - I) 



(3.10) 


y = [>£( 0, 0). y F ( 1,0) y £ (* - 1,0)1 = 


(3.11) 


3.1.6 


unstack 


Format : 

[w^unstackiy, stppbtp > 

Description : 

Convens an ETIS vector y into a traditional discrete 

Inputs : 

y An ETIS vector of the form of Eqn. (3. 1 1 ) 


stppbtp See Table 3.1 

Outputs : 

w A vector of the values of y(m,n) as in Eqn. (3. 10) 


vector w. 


NOTE: Additional M-Files used 


in support of those described in this section are described 


in Appendix B. 
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3.2 SIMULATION M-FILES 


3.2. 1 mrsim 

Format: 
Description : 


ly, 1 1 =mrsim(plt, nplt, cmp, ncmp, su. sy , sz, stppbtp , simpstp. inputl . outputl. ur. stp, XO) 

Generates a ume response simulation of a continuous or discrete plant in feedback with a mult, rate 

nZtnT f r ° UUne 25511,1,65 P0SU,Ve feedbaCk “ Sh0W " ^ W,th thc ‘on'-nuous inputs 
nd outputs of interest sampled or updated every stp/simpstp seconds. This configuration allows 

the user to simulate the inter-sample response of the multirate system. 



Tf = stp/simpstp 


Inputs : 


Outputs : 


pit, nplf A description of the plant, see Table 3.1. 

cmp. ncmp. su. sy, sz, stppbtp .A description of the mult, rate compensator, see Table 3 I 
simpstp This routine can calculate the response of the closed-loop system between the 
sample/update activities when the plant ,s continuous (stp > 0). simpstp is the number of 

:7::t of ,he — simpstp - - - — -- ^ * an d 

inputl A vector specifying which plan, inputs are connected to the compensator outputs For 

example mputl =[14] indicates that the compensator outputs are connected to the firs, and 
fourth plant inputs. 

outpu.1 A vector specifying which flan, output, are connected to the compensator inputs Fo, 

example oulputl = 12 4| tnd, cates that the compensator tnputs are connected to the second and 
lourth plant outputs. 

* “r T "7 COn ' a '" , " g ,h ' Va '“' S ° f mp “' “ evtry "me fo, 

K stppbtp simpstp sample times, ur has the form 

«• = I tt^tO.O). iyO,l) Up (k - 1 , stppbtp simpstp - I)] (3 )2) 

Note that in Eqn. (3.12) the sample/update period is one stp/simpstp. This sample period should 
not be con used w„h the sample period of the multirate compensator which is simply one stp. 

The second index ,n the double index notat.on in Eqn. ,3.12) takes on values between zero ard 

simpstp- stppbtp . 

stp see Table 3.1. 

XO A vector containing the initial conditions of the plant state vector. If XO is omitted, the initial 
conditions are assumed to be zero. 

y A vector containing the response of the plant every ( stp/simpstp ) time, y has the form 
y = I Vp(O.O). Vp,0. 1 ). ,\p(k-l. Stppbtp simpstp - I )| 


(3.13) 
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Note that the second index in the double index notation corresponds to the second index of 
Eqn. (3.12) 

t a vector containing the times corresponding to the columns of y. When stp = 0 then t indicates the 
number of sampling periods. 

Comments: This M-File uses the ETIS of the system to compute the time domain response. Hence, if k in 

Eqn. (3.12) is not a multiple of simpstp-<V where </V=BTP/STP) the system response is only 
simulated to the nearest multiple of simpstp-zV 

3.2.2 mrfeedback 

Format: (ae. be, <*. de ) =mrfeedback( pit, nplt, cmp, ncmp, su. sy , sz.stppbtp. input!, output! , stp) 

Description: Creates a closed-loop ETIS system from a discrete or continuous plant, depending on stp, and a 

multirate compensator using positive feedback. 

Inputs: pit. nplt The plant description. See Table 3. 1 . 

cmp. ncmp, su. sy , sz. stppbtp The multirate compensator description See Table 3.1 . 
inputl A vector specifying which plant inputs are connected to the compensator outputs. For 
example inputl = [1 4) indicates that the compensator outputs are connected to the first and 
fourth plant inputs. 

outputl A vector specifying which plant outputs are connected to the compensator inputs. For 
example outputl = [2 4] indicates that the compensator inputs are connected to the second and 
fourth plant outputs 

stp See Table 3.1. If stp is omitted mrfeedback assumes stp=0. 

ae, be, ce, de The system matrices o i the ETIS closed-loop system with the form of Eqn. (3.6). The 
inputs and outputs of this system are the inputs and outputs of the discrete plant 


Outputs : 
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33 SYNTHESIS M-FILES 

The synthesis algorithm is implemented as a series of script files. They allow the designer to optimize the 
performance of the closed-loop system shown in Figure 3. 1 by calculating values of A z , B z% C z> and D z , such that the 
cost function of Eqn. (2.6) is minimized. There are three distinct steps to the optinrzation process: 1 ) entering the 
data which describes the problem; 2) discretizing and preprocessing the data; and 3) optimizing the compensator 
gain values. Each step of the optimization process is described in detail in the following paragraphs. 

Although the synthesis algorithm was designed to solve the multirate sampled-data problem, it is fairly general 
and can be applied to a variety of related problems as well. It can compute optimum compensator gains for a 
sampled-data system consisting ot a continuous plant in feedback with a discrete compensator as shown in 
Figure 2. 1, or for a discrete system consisting of a discrete plant in feedback with a discrete compensator. In either 
case, the compensators may be single-rate or multirate and may have either time-invanant or periodically time- 
varying digital processor gains. 

33. 1 Input V ariabies 

The synthesis algorithm looks for specific Matlab “workspace" variables to define the optimization problem. 
(A workspace variable is defined from the » prompt as opposed to being passed as an argument to a function.) 
These variables define things such as the plant dynamics, the cost function and the compensator structure and are 
defined in Table 3 2. The user must assign values to these variables before beginning the optimization. This can 
either be done manually or automatically from a script. An example of a script which defines these variables is 
provided in Appendix D. 




Plant 

Continuous or Discrete 


Vc 


/s 


7 0H 


u(mn) 


Zero Order Hold(s) 


"itmrHl ) - A^ntfrrxn) ♦ B^yimn) 
= C/riftmri) + D/n)^rr\n) 


%mn) 



Sampler(s) 


*rr\n) 


Digital Processor 

Figure 3 t Closed-loop MuUirute System 




Np 

pit 


wxx 


^cVon^T° fPlamPertUrbati0nS NP= ‘ indicateson| y ,hen °m>nalpant,s defined. Sec 
State-space matnces of the plant(s). Given the I th continuous plant perturbation 

|>„(0j H ' X ' {!) + J ‘\ o | (3.14b) 

or the « * discrete plant perturbation, depending on the value of stp defined later. 

x f (m,n+l)= F i x i (m.n)+ G, u i (m.n)+ Hy(m.n) (3.J4C) 

ty«(«.*)J + o f (3.14d) 

vvhere u is the control input; w, is zero mean white no.se input.; y a ,s the enter, on output; 
and y„ are the outputs which are sampled by the compensator. Then 


Pit: 


rr p , 

Ft Gj 


F Np Gn p T 

i — 

a: 

h 2 J 2 


_ W Np Alp 


(3.15) 


For convenience we assumed in Rqns. (3.14) that the input vector ,s part.t.oned into 
control input u, and the disturbance input w,. The software, however, can accommodate 
a more general form for the input vector. Using the variables ucol and ncol. defined later 
t e user can specify which inputs correspond to control and noise inputs respectively. Anv 
row o t e input vector can be either control, disturbance or both. (We interpret an input 
vector specified as both control and disturbance as two separate inputs, one a control input 
and the other a disturbance input whose distribution vectors have the same values , 

Similarly, the outputs >■„ and are specified by crow and srow respectively, and 
any row ot the output vector can be criterion output, sensor output or both 
NOTE: For a discrete plant, the process noise distribution matrix is unity. 

A vector ot the PSD values ol the white noise disturbance n / If 

W i 6 (T)=E[w,{t- nwf(/,) 13.16.) 

<,nf Wj(6(k) + 6(j)) = E{K,(ni.n)*l(m + k.n + j)) <3. 16b) 

depending on the value of stp, then 



continued on follow in# page... 


__J^li£^^^ rks P ace Variables Defimng Optimization Problem (continued from pre nous page , 
Variable Description 

vr The PSD values of the discrete sensor noise covariance. If the output, sampled every STP. 

is 


where 

then 


>«<* ) = H,x,(k) + vpk) (3.17) 

Vj8(j)= E{vi(k)v[(k + ;)} 

vt=[V\ V 2 ... V Np ) 


A vector specifying which columns of G, correspond to w,. If ncol = ( 1 4). then the plant 
disturbance vector w is comprised of the first and fourth plant inputs. 

A vector specifying which columns of G, correspond to u,. If ucol = [2 4], then the plant 
control input vector u is compnsed of the second and fourth plant inputs. (Again, we 
interpret an input vector specified as both control and disturbance as two separate inputs. 

one a control input and the other a disturbance input whose distribution vectors have the 
same values.) 


srow A vector specifying which rows of C, correspond to y st . If srow = ( 1 3], then the first and 

third continuous plant outputs are sampled and connected to the compensator. 

crow A vector specifying which rows of C, correspond to y ci . If srow = [13], then the first and 

third continuous plant outputs are used to calculate the cost given in Eqn. (3.18). 

qla, q2a, maThe cost function weighting matrices for either a continuous or discrete cost function, 
depending on the value of stp defined later. The cost function has the form 


^ continuous = 

t — >oo 


1 = 1 


Jf.T/,, T,..] 

'01, 

Mi' 

yd | 

|[3ci'(0 li i (t)J 

mJ 

02i. 



(3.18a) 


or 


■/ discrete ” ^ irn 
m— »oo 


Stppbtp 


Slp^SHp- 

1=1 11=0 


| y Ci (m,n) 

T 

01 1 M, Tv t .,(m.u)]| 



L Mj 02,J|_ u,(m.n) | 


(3.18b) 


where i corresponds to a particular plant condition 


qla = ( (2l I 012 Gl Np] 
q2a =(02i 022 ... 02Npl 
ma = ]/V| N 2 ... Visjp] 


qa and ra must be symmetric and positive semi-definite. 
If nta is left undefined its value is assumed to be zero. 


continued on followtn r page . 
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Table 3. 

2 Mailab Workspace Vanables Defining Optimization Problem ( continued from previous page ) 

Variable 

Description 

emp 

The compensator gain matrices as defined in Table 3.1. The digital processor gains can be 
either periodically time-varying or time-invariant. 

IMPORTANT: The user must provide an initial guess for emp which stabilizes ail plant 
perturbations. 

empfree 

This matrix specifies which of the gains in emp are fixed, empfree has the same 
dimensions as emp, but its elements are either 0 or 1. AO indicates that the corresponding 
element in emp is fixed and should not be optimized. A 1 indicates that the corresponding 
element in emp is free and can be optimized. 

su, sy , sz 

The sampling schedule description. See Table 3.1. By selecting appropriate values for su. 

sy, sz and Stppbtp the designer can specify whether the compensator is multirate or single- 
rate. 

Stppbtp 

The same as stppbtp in Table 3.1. Note that the first letter of Stppbtp is upper case, 
indicating this is a global workspace variable. See Section 3.3.2. 

stp 

STP as defined in Tabic 3. 1 . 

If stp > 0 the algorithm assumes that the plant (pit), the processor noise (wxx) and the 
cost weighting matrices (qla. q2a. and ma ) are all defined for a continuous 
plant in feedback with the discrete compensator. 

If stP = 0 the algorithm assumes that the plant (pit), the processor noise (wxxi and the 
cost weighting matrices (qla, q2a. and mat are ail defined for a discrete plant 
in feedback with the discrete compensator, and that the sampling period of the 
discrete plant is one STP of the multirate compensator 


3.3.2 mroptjnit 


Format : 


mroptjnit 


Description : 


Inputs : 
Outputs : 
Commei.r. 


This script operates on the workspace variables defined in Table 3.2. It discretizes (if necessary) the 
plant, process noise and cost function; and defines a set of global workspace variables used by the 
optimization script ( mroptoptim \. This script must be rerun after making an) changes in the 
workspace variables defined in Table 3.2. with exception to changes in the values of emp or 

empfree . 

mroptjnit is a script which uses the workspace variables defined in Table 3.2. 

mropt_init save it's output to the global workspace variables defined in Appendix A 

The global workspace variables defined by mroptjnit are used by both the optimization scr.pt 
mropt.optim and by the M-File calc LQGcost All global variables generated by mroptjnit 
begin with a capital letter. This helps to differentiate them from other workspace variables. A 
brief description ol these global variables is given in Appendix A. 
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3.3.3 mropt_optim 

Format : 

mropt_optim 

Description: 

mropt_optim computes optimum values of the d.g.tal processor gains such that the cost function 
defined in Eqn. (3.18) is minimized. 

Inputs : 

mropt_optim uses the global workspace variables generated by mropt Jnit in addition to the three 
workspace variables uefined in Table 3.3. 

Outputs : 

The script returns a value of the cost function, the gradient of the cost function with respect to the 
tree compensator parameters, and the optimum values of the digital processor gains. These 
values are '.utomatically displayed on the computer screen and simultaneously stored in the 
workspace variables defined in Table 3.4. 

Comments: 

1 ) As already noted, if the user modifies the problem definition by changing the values of the 
variables defined in Table 3.2, the senpt mropt Jnit must be rerun before running nropt_optim 
The exceptions are changes to the variables defined in Table 3.3. 

2) mropt_optim can be computationally intensive. The user can abort the optimization with a 
Control -C key sequence. Also see mropt.extract. Section 3.3.4 


Variable 


Table 3.3 Input Workspace Variables for mropt.optim 
Description 


This variable contains the compensator gains and is the same variable as described m 

Table 3.2. These are the starting values for the optimization and must stabilize all plant 
perturbations. 


cmpfree 


This matrix specifies which of the gains in cmp are fixed and is the same variable as 
described in Table 3.2. 


During optimization the compensator gains are scaled using a variable named dscale to help 
improve the numerical accuracy of the search. They are scaled as follows 

Xscaled - inv(dscaJe) Xunscaled 

where Xunscaled is an unsealed vector containing the free compensator gains and Xscaled 
are those values scaled. Typically dscale is a diagonal matrix chosen so that all the 
elements in Xscaled have the same magnitude. The user can either manually set dscale or 
have mropt_optim set it automatically. If dscale exists and has the correct dimensions 
then mropt_optim uses that value of dscale. otherwise a new value for dscale. based on the 
current values in cmp, is automatically calculated. Typing "claar d.cal* before 

running mropt.optim forces the routine to scale the problem based on current values in 
cmp. 
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Table 3.4 Output Workspace Variables tor mropt_optim 


Variable 


Description 


cmpf 


jcost 


djdp 


This contains the optimized values of cmp after 
optimization. 


This contains the value of the cost function after 
optimization. 


cmp after mropt_optim has completed the 
mropt_optim has completed the 


This contain . i scaled value of the gradient of the cost function with respect to the 
compensator values after mroptoptim has completed the optimization. 


3.3.4 


mropt_extract 


Format : 
Description: 

Inputs : 
Outputs : 

Comments : 


mropt_extract 

Alter aborting the script mropt_optim. usually w.jh a Control-C key sequence, the script 
mropt_extract can extract the value of the “optimized” compensator at the last iteration. 

None. Can only be rerun after aborting mropt_optim. 

The optimized values of the digital processor gains at the last iteration before mropt_optim was 
aborted are displayed on the computer and stored in the workspace variable cmpf. 

Sometimes it is necessary to abort mroptoptim before is has completed -he optimization. The user 
might wish to rescale the free parameters by changing dscale. or might decide to free up 
additional gains in cmp by changing empfree. The user can restart the optimization from the last 
iteration by setting cmp=cmpf and rerunning mmpt_optim. 


3.3.5 
Format : 
Description 

Inputs : 


Outputs : 
Comments: 


calc„LQGcost 

(jcost h= calc_LQGcost ( j ) 

Calculates the discrete LQG cost assuming a sampling period of one stp. This roufne uses the 
global workspace variables computed by mroptjnit Use th.s routine to find the minimum value 
ot the cost function attainable by mropt_optim 

j = 0 or is undefined, then the routine calculate the sum of the LQG costs for all Np plant conditions 
corresponding to the cost J discrete given in Eqn. (3.18) 

0<j<Np then the routine calculates the LQG cost for only the j* plant condition 
jcost The value of the LQG cost function 

calc_LQGcost uses the Controls.Toolbox routines dlqr and dlqe. They require that ra in Table 3.2 
be positive-definite. 


NOTE: Additional M-Files used in support of those described in th.s section arc described in Appendix B 
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Section four 


Synthesis Example 


4.0 PROBLEM DESCRIPTION 

In this section we present an example design problem which utilizes the M-Files discussed in the previous 
section. Our objective is to design a multirate compensator for a lightly damped mass-spring-mass system. 

The mass-spnng-mass (MKM) system is shown in Figure 4.1. It consists of two masses connected by a spring 
and damper. The control inputs are the force inputs and \ the sensor outputs are the displacements x \ and 
and the disturbance input is w. The sensor outputs are corrupted by discrete sensor noise (not shown on figure) with 
covariance v. Nominal values for the plant parameters, along with those for one known plant perturbation, are given 
in Table 4.1. These parameter values result in a system with two poles at the Oiigin associated with x\, and two 
lightly damped poles associated with X2 The open-loop poles are listed in Table 4.1. A M-File which defines the 
system matrices for the MKM system is given in Appendix C. 

The goal of our design is to increase the damping coefficient on the lightly damped poies to £=() 707 without 
shifting their natural frequency, and to move the two poles at the origin so that the x\ response has a frequency of 
0.825 rad/sec with a damping coefficient C=0.707. 




/ 
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4V/ //////.'■ /// 


Figure 4.1 Mass-Spnng-Mass System 
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Table 4 1 Mass- Spring- Mass Parameter Values 



Nominal Plant 

Perturbed Plant 

M 

! 0 kg 

! .0 kg 

m 

0.1 kg 

0.2 kg 

k 

0.01 N/m 

0.01 N/m 

b 

0.01 Ns/m 

0.01 Ns/m 

E( w T (f-t)vi(n( 

O.l&riN 2 

o.i&rtN 2 

E{ v^ik+j )iU)} 

0.001 &/)m 2 

0.001 &/>m 2 

Open-Loop Poles 

0. 0. -0.055± 3.316i 

0. 0. -0 030t2.449» 




4.1 COMPENS/. TOR DESIGN 

A " ~ ^ - « * < - *. ^ 

computer, instructional comments are in italic. * ^ Plain Cexc ,s returned by the 

S " P ‘‘ , '“ d,a “ C ' 0, ‘ ta *■>»• <M> »hich characterize 

we have two plants - T nomml and^he ^ ^ Piant pertUrbatioa In our c ** 

function by designing two continuous time LQ regulators ^£*77*,** <V 0Ur c ° s ' 

petturbed plants in the desired locations. The coa functions arc lh « ulosed-loop poles for the notntnal and 


^nominal = E{5.5 jt, 2 + 2.2x\ + lOOOOw 2 + 10« 2 } / 10.6 
'perturbed = E{6-5* 2 + 4.8i 2 + lCXOu 2 +10u 2 }/ 8.0 


(4.1a) 
(4.1b) 

where the costs have been scaled so that the mp . r 

performance of our multirate compensator ,s th^sum of the th “ ^ C ° $t furc,ion characterizing the 

'perturbed) the optimization routine will attempt to find a single 'mmae ^ fU " CU ° nS ' By minimi ’ in g (^nominal + 
the LQR designs for both the nominal and the i^rt h w P^sator w ose performance approaches that of 

exist. Instead the resulting compensator will represent a P ^ ^ ^ ^ $UCh 3 com P*nsator to 

Plant condition and good performance at the pert^d cond 2^ ^ Perf ° rmanCe 31 — 1 

^ 2: ^Z r I ri “' C0 " PeM * ,Or *" d — * *— - .he desired eh«d- 

- -rcssir;:,::: :^ h “ m,h * and - ^ — w,u,z 2 

order loops, one from X! to u\ and another n to in aT. 3 Th ^ SlnJCtUre - If cons,s:s of tw « coupled first- 
x| ton, loop is sampled and updated every 0~8 seconds whT hT™™ ° f ,hlS SUucture ,s shown In Fi ? u re 4.2. The 
often, every 0.2 seconds. These sampling rates ar 1 e 1 - x 2 1° U 2 loop is sampled and updated lour time as 

loop. In addition, the compensator accounts forTTeTTf 3 ^ T closecHoo P bandwith for each 

schedule is shown in Figure 4.3. * f Sampl ‘ ng penod computational delay. It’s sampling 



Figure 4.2 Block Diagram of Closed-Loop Multirate System 


*1 Command 
*2 Command 
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STP=0,1 sec — 


^■+111111? Sampler 

I t I Y I ? Sampler 

H4 - + I I I I II ZOH 

- I- ? I ? I ? I ? I ZOH 

^siow h - t-4 I I I I M State Update 

^ast ?-+ ? I ? I ? I ? State Update 

c — J L l 

BTP=0.8 sec 


*1 

*2 

u i 

u 2 


U 


? = sample/update activity 
Figure 4.3 Sampling Schedule for Muitirate Compensator 


Step 3: Dedg. . c.mp. 0 s. Wr »,«, ,h« sunplto, sch^dul, .„ d s ,„ cm „ s „ p 2 (h> , - 

plant perturbations. This will be the starting point for the optimization. 

We designed a compensator using successive loop closures and root locus for the nomtnal plant The 
compensator has the same structure and sampltng schedule as the multirate des.gn discussed in Step 2 except there ts 
no input coupling and we did not account tor computational de.ay. This compensator does ho^„« 

^ Pl3ntS eVe " WUh ^ ~ “* S ° « - « - starting “die 

Refer to Reference 1 for an introduction to successive loop closures. 

Step ™ .3 4 LOad the pr ° b,ein definition int0 MaUab’s workspace variables defined in Table 3.2 

We defined the workspace vanables using the Matlab scnpt provided in Appendix D. To load the workspace 
variables for this example, type woricspace 

»mropt_ 


Step 5: Initialize the workspace variables and generate the 

Type 

»mropt_ini t see Section 3.3.2 


necessary global variable using mroptjnit 


At this point you can compuf he minimum value of the cost function as follows (see Section 3.3.5) 

»calc_LQ<3cont see Section 3.3.5 
» ans = 

2 . 0 1 0 8 e + 0 0 this is the minimum cost our design could achieve 
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Step 6: Calculate optimum values of the digital processor gains using mropt.optim 

A partial output of the optimization follows 


»aropt_opti» see Section 3.J.3 


y selecting dscale a diagonal dscale is automatically selected if the variable 

dscale does not exist, or is not compatible with the 
1 current compensator 

a culate gradient only? (y or n)n if we replied y the algorithm would compute only 

the gradient of the cost with respect to the free 
digital processor gains and the value of the cost 
for the current compensator gains in emp 


f- COUNT = 
FUNCTION = 


the total number of function evaluations 

'dZruZ!Z£!! e costfunction at the current deration. Refer to Matlab FMINU 


f -COUNT 
2 

-> Plant 
-> Plant 
-> Plant 
-> Plant 
-> Plant 
-> Plant 
-> Plant 
-> Plant 
-9.9655e 


FUNCTION 

176 

unstable 
unstable 
unstable 
unstable 
unstable 
unstable 
unstable 
unstable 


-01 


STEP-SIZE GRAD/SD 

0.01 -9 . 48e+07 

while calculating J 
while calculating J 
while calculating J 
while calculating J 
while calculating J 
while calculating J 
while calculating J 
while calculating J 


LINE-SEARCH 


The algorithm encountered a destabilizing 
solution. It automatically adjusts the step 
size and continues the optimization 


-1 . 63 99e-02 


1 . 0030e+00 
-1 . 0676e-02 
-2 . 2977e-02 
-9 . 4363e-01 
-1 . 3703e-03 


these are the values of the scaled digital processor gains 
at the current iteration 


1 . 0529e+00 
1 . 0030e+00 
1 . 0036e+00 
1 . 0529e+00 
1 . 0085e+00 


? . . 14.2 2.45e-i3 int_st 

irni zat ion Terminated Successfully 
Gradient less than options (2) 

NO OF ITERATIONS=314 
-8 . 0984e-01 
- 1 . 6700e+01 
1 . 8697e+00 
2 . 2849e+01 
7 . 4219e-01 
-7 . 3479e+00 
-8 . 6882e+00 
7 . 9308e+00 
1 . 0706e+00 
3 . 8424e-03 
1 . 0367e+00 
1 . 4037e+00 


Final Results 


Gain 

Gradient 

-1 . 6197e-01 

-4 . 6527e-06 

-1 . 5700e+01 

6 . 6547e-04 

1 . 8697e+00 

-2 . 8675e-04 

2 . 2 849e+01 

4 . 5050e-04 

7 . 4219e-01 

1 . 6 57 6e- 04 

-5 . 5109e+00 

-8 . 9955e-04 

-8 . 6882e+00 

-4 . 8353e-05 

7 . 9308e+00 

-8 . 5171e-04 

8 . 5648e-02 

-1 . 0838e-04 

1 . 9212e-03 

6 . 9640e- 04 

6 . 2202e-01 

3 . 4133e-03 

1 . 4037e-01 

3 . 2564e-04 


Final cost - 4.62482 the optimum value of the cost function 

Optimized compensator gains 
empf = 


1 . 9212e-03 
0 

8 . 5648e-02 
0 


0 

1 . 4037e-01 

0 

6 . 2202e-01 


1 . 8697e+00 
2 . 2849e+01 
-1.6197e-01 
-1 . 6700e+01 


-8 . 6882e+00 
7 . 9308e+00 
" • 4219e-01 
-5 . 5109e+00 


Elapsed time: 1516.05 sec 
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4.2 DESIGN ANALYSIS 

We looked at two criteria when evaluating our mullirate design: 1 ) the tinal cost; and 2) the step response to a 
command input. 

4.2.1 Final Cost 

The final cost for the multirate compensator is calculated automatically by the optimization routine 
mropt.optim . The LQG cost is computed by the function calcJLQGcost in step 5. For our system 

^multirate = 46 and -/LQR = 2.0 


4.2.2 Step Response 

We obtained the response to a unit command step to jc| as follows 
» [a, b, c, d] anricm; get the nominal plant 

»b» [b zeros (4,1) ] >c«[c; [1 0 0 0]]; create a reference input 
»d«[d zaroz(4,l);0 0 -1]; 

»ur« [ zeros (2,100); ones (1,100)]; generate a step input function for the input 

and compute the time domain response, see Section 3.2.1 
»Iy,t]-mr«im( [a b;c d] , 4 , cmpf , 2, au, ay, az, Stppbtp, 1, [1 2], [5 3],ur,atp); 

assuming cmpf. su.sy.sz. Stppbtp, and stp were previously defined and 
that cmpf contains the optimized values of the digital processor gains 

»plot(t,y( [1 3] , * ) ' 

We similarly obtained the step response for the LQR design. The results are shown in Figure 4.4. 

As expected the performance of the multirate design does not match that of the LQ regulator. The LQ regulator 
was optimized for only the nominal plant. The multirate designs on the other hand represent a compromise, 
stabilizing both plants and providing “optimum" performance for both. In addition, the LQ design was a continuous 
full-state design whereas the multirate design was discrete second order and used only two state measurements. 
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appendix a 


Global Workspace Variables 

the first letter of every global 
For each plan, perturbation. the discrete closeddoop system can be written tn the form: 


The follow, ng global workspace variables are defined by mropl inrt Note that 
vanable is capitalized. 


*i(/n,n+l) = FjMmM)+ G t u(mM)+ Ww t 
V|(m,/i)= Hjx(mM)+ ttv, 

u,(m.n)= [51(«)P( W )52(> 1 )+ S3(n)]y s (mji) 


where 


F,= 


r F t 0 
0 0 


Gj = 


G, 0 
0 / 


. w = 


r / oi 

0 0 


W,= 


H, 0 

0 / 


l/ = 


0 / 
0 0 


, P = 


D, c. 
B . A - 


(A.l) 

(A.2) 

(A.3) 

(A.4) 


and w is a discrete zero mean white noise input with 

E{ w, < m, n )wj( m + j , n + k )} = [S(k ) + S( j )} R, 


(A.5) 


tC fd' T"°" *" d ,h ' 0Verbm ' d mamc “ " Pi- matrices. The 

-me. — — - — — 

The discrete cost function tor the system in Eqns. (A.l -A.3) can be 


written as 


Np A/-I f 


.t,(m.rt)l' 

Qli Mi lf.r,(m.n)l| 


.M, (kiUm m.n) I 


(A.6) 


See Reference 1 for more information. 


The global Matlab workspace variables are then defined as follows. 


Plant Description: 

Np = number of plant conditions 

MF = [F\, A : A Np | NcF = 

MG= [ G \ , G 2 O'Np | NcG= 

( W| , Wi IV'Npl NcWs 

MH=[//|,// 2 // Np j NcH= 

V'Np] NcV 

MR* (/?, . /?2 /? Np | NcR 


number of columns in h, 
number of columns in G, 

• number of columns in VP, 
number of columns in //, 
number of columns in V, 

■ number of columns in R , 
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r ^mtm n tk 


t+Ql Ft-WXj 


Compensator Description: 

Nzhcmp = number of states in the periodically time- varying model of the multirate compensator 
Stppbtp = the global version of stppbtp 

Slk = (51(0), 51(0) 51(AM)] Nsl = the number of columns in 51 

S2k = ( 52(0), 52(0) 52W- 1 )] Ns2 = the number of columns in 52 

S3k = (53(0), 53(0) 53(AM)| Ns3 = the number of columns in 53 

P*v: abs( Ptv> = number of columns in the compensator gain matrix [A z . B z , C,, D z \. Ptv is positive if the 
digital processor gains are time-invariant and negative if they are periodically time-varying. 

Cost Function Description: 

MQ1- [ Q\ j, ~ Q\ N p ] NcQl= number of columns in Qn 

MQ2- [Q 2 \, Q22> ••• Ql Np) NcQ2= number of columns in Qi t 

MM= [M i , M 2 , ... Af Np J NcM= number of columns in M, 

Miscellaneous Definitions: 

Ppn , Ppm, and Ppinit : used internally 

Ppfree = index to the '‘tree" digital processor gains 

Fd - closed-loop state transitions matrix of the last computed ETIS system 

Xc_global = value of the compensator gains at the last iteration 

Unstable . If Unstable =0 the compensator gains at the current iteration stabilize the plants 
Otherwise they destabilize the plant 


appendix B 


Support M-Files 


B.l MODELING SUPPORT M-FILES 


B.I.l mrspiit 


Format. 


[a, b, r, d ] = mrsplit(P, n ) 


Description : Separates the composite gain matrix P into its constituents, where 

„ I* cl 


Inputs : 


Output 


Psee Eqn. (B.l) 
n the number of columns in a 

a, b,c, d the gain matrices in Eqn. (B.l) 


mnmakesk 


Format : [slk, nslk, s2k, ns2k, s3k, ns3k]=mnnaker w su, sy, sz, stppbtp ) 

Description: The periodically time-varying representation of a multirate compensator in Eqn. (2.2) can be written 


y(m.n) ]_[/>(«) C(n 
x(m.n + 1) B(n) A(n 


)Tn(m,n)] j D.{n) C-(/i)l i["w(m.n) 

= \ Sl(n) - - S2(n) + 53(/i)> 

)J[jc(m,n)J ^ B.(n) /W/i)J j[jc(m.n) 


where A : , B z . C- and D- are the digital processor gains and 51. 52 and 53 are switching matrices 
which model the compensator’s sampling schedule, mrmakesk creates the periodically time- 
varying switching matrices given the sampling schedule description. 

Inputs : su. sy , sz . stppbtp see Table 3. 1 

Outputs : slk = [51(0) 51(1) ... 51(stppbtp-l )[ 

s2k = [52(0) 52(1) ... 52( stppbtp- 1)| 
s3k = [ 53(0) 53( 1 ) ... 53( stppbtp- 1 )| 

nslk. ns2k and ns3k the number of columns in 51( ), 52( ) and 53( ) respectively 


Format : 


mrgetsk 


(sl.s2. s3)=mrgetsk(k. slk. nslk. s2k. ns2k. s3k. ns3k) 


Description: Extracts the individual switching matrices for a given STP from slk. s2k and s3k 

Inputs : k an integer specifying for which STP the switching matrices are to be extracted 

slk. nslk, s2k. ns2k. s3k. ns3k see Section B 1 .2. 

Outputs: si. s2. s3 switching matrices corresponding to 51(k). 52(k) and 53(k) respectively. See 

Section B.1.2. 
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B.2 SYNTHESIS SUPPORT M-FILES 


B.2.1 calc_djdp 

Format : djdp= caic_djdp( p guess, dscale ) 

Description: Used by the routine mropt_optim to calculate me derivative of the cost J with respect to the free 


Inputs : 


Comments: 


compensator parameters. See Section 2.3. 

pguess a vector of the values o t the free digital processor gains at the current iteration 
dscale see Table 3.3. 

This routine uses the global workspace variables generated by mropt_init 


B.2.2 calcj 


Format : 
Description : 
Inputs : 

Comments : 


djdp= calcj ( pguess, dscale) 

Used by the routine mropt.optim to calculate the cost J at the current iteration 

pguess a vector ot the values ot the free digital processor gains at the current iteration 
dscale see Table 3.3. 

This routine uses the global workspace variables generated by mroptjnit 


B.23 calc_L 


Format : 
Description : 

Inputs : 


[lacc]= calc_L(pk, f, g, w, h, ql, q2, m) 

Called by calc_djdp to calculate the steady-state values of a Lagrange multiplier. (The multiplier is 
used to adjoin the stability equality constraints to the cost function.) 

pk gams for the periodically time-varying representation of the multirate compensator 


pk = ( P(0), P(\ ) P{ N - 1 )] and P(i) = 


DU) 

B(i) 


C(i)' 

Mi) 


(B.3) 


0(0, C(i), B(i), and Mi) are the state space matrices for the periodically time varying 
representation of the compensator. See Section 2.2 and Section B. 1 .2 
t g, w, h plant description matrices for the current plant perturbation corresponding to F,, G„ W and 
Hi in Eqn. (A.4). 

ql. q2. m cost weighting matrices for the current plant perturbation corresponding to (?i„ Q 2 , and 
Mi in Eqn. (A.6). 


Output: 


Comments: 


lacc periodically time varying lagrange multipliers where 
lacc = ( A(l) A(2) ... A(N-\) A(0)] 
and A(i) is the lagrange multiplier for the i lh STP 
This routine uses the global workspace variables generated by mropt.init 


B.2.4 


calc_P 
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Format : 
Description : 




[pk]=calc__P('p) 

Called by calc_j and calc.djdp to compute the periodically time-varying representation of the 
multirate compensator 

p the digital processor gains in the form 


p = (p(0),p(l) p{ A/- 1 )] where p{i) 


D z (i) CM) 
B z (i) A z (i) 


Outputs : pk matrix of the compensator gains for the periodically time-varying representation of the multirate 

compensator 

pk = (P(l), P(2), PO),... P(A0) where °(i)= Sl(i)p(i)S2(i)+ S3(i) 

See Eqn. (A. 2) or Section B.1.2. 

Comments: This routine uses the global workspace variables generated by mroptinit 

B.2.5 calc_X 


Format: 


[xacc |= calc_X(pk, f, g. w, h. v. r i 


Description: Called by calc.djdp and calcj to calculate the steady-state covariance values of the closed-loop 

system 

Inputs : pk matrix of the compensator gains for the periodically time-varying representation of the multirate 

compensator, see Section B.2.4. 

t g, w. h plant description matrices for the current plant perturbation corresponding to F t . G,. W. H, 
and V in Eqn. (A. 4) respectively 

r process and sampling noise covariances for the current plant perturbation corresponding to R, in 
Eqn. (A.5). 

Comments: This routine uses the global workspace variables generated bv mropt init 

B.2.6 disc_cost 


Format : 
Description: 


[qld, q2d, md)=disc_cost(plt. nplt, stp, ql, q2, m) 


Computes the discrete cost function weighting matrices such that the cost associated with a 
continuous plant in feedback with a single-rate compensator is identical to the cost associated 
with the discretized plant in feedback with the single-rate compensator. The continuous cost 
function is given in Eqn. (2.6); the discrete cost function is given by 


J di 


iscrete 


= lim 


f 1 

r roid 


^mm.ni 

_A/J 



(B.4) 


disc.cost assumes the control input is updated using a zero-order-hold. 

Inputs: pit, nplt stp see Table 3.1 

ql , q2. m the continuous cost function weighting matrccs in Eqn. (2.6) 

qld, q2d, md the discrete weighting matrices in Eqn. (B.4) 


Outputs : 
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B.2.7 dtsc_noise 

bormat : [w ^disc_noise(wxx, pit, nplt stp ) 

Description. Computes the covariance of a discrete-time white noise process disturbance such that the expected 
value ot the states of a continuous plant, and the expected vJue of the states of the discretized 
plant are equivalent at the sampling instances. Thus, given the continuous plant 

*(/) = A.r(f) + G H .w(r), where E(w(t)w T U + r») - 5(r)wxx 
and the plant discretized with a zero-order-hold at a sampling period T 
x d( m. n + 1 ) = Adx d ( m, n ) + w d ( m, n ), where £( w d ( m, n )vt j(m + k . n + j ) ) = < S( k ) + 8< j ))w 
then E{x d (m,n)x])(ni,n))= E{x((mN + n)T)x T ((mN + n)T)} 

Inputs : wxx PSD of the continuous processor noise 

pit, nplt, stp see Table 3. 1 The inputs to pit are assumed to be only the process noise 

Outputs. w PSD ot the discrete process noise 


B.2.8 diyap2 

Forhiat : x= dlyap2 ( a, c, method ) 

Description : Solves the discrete lyapunov equation 

x= axa T + c 


Inputs : 


Outputs : 


a, 

if 


c system matrices in Eqn. (B.5) 
method = ‘Bartels’ then the algorithm solves 
otherwise Eqn. (B.5) is solved 


Eqn. (B.5) using Bartels method 
using partial sums 


x the solution to the lyapunov equation 


(B.5) 


B.2.9 get_cost 


Format : [ql , q2, m ^=get_cost(i ) 

Description. Extracts the cost function weighting matrices Ql, Q2 and A/ for the i ,h plant perturbation from the 
global workspace variables MQI, MQ2, and MM See Appendix A. 

Inputs . i an index to the desired plant perturbation 

Outputs : ql , q 2, m the cost weighting matrices corresponding to 01 , 02 and M respecuveiy 

Comments. This routine uses the global workspace variables generated by mropt Jnit 
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B.2.10 get 

Format : 
Description : 

Inputs : 

Outputs : 

Comments : 


_plt 

[f, g, w, h, v,rj=get_plt(i) 

Extracts the plant matrices F„ G„ W, H„ V. and Rj for the i ,h plant perturbation from the global 
workspace variables MF, MG. MW, MH, MV and MR. See Appendix A. 

i an index to the desired plant perturbation 

f, g, w, h, v, r plant description matrices tor the plant perturbation corresponding to F t , G t> W, //,*, 
V % and /? f respectively 

This routine uses the global workspace variables generated by mropt.init 


B.L11 

Format : 
Description. 
Inputs : 
Outputs : 


get_ppfree 

ppfree = get_ jjpfree ( cmpfree , ncmp ) 

Called by mropt_optim to determine which compensator gains will be optimised 
cmpfree. ncmp see Table 3.1. 

ppfree vector whose elements point to the free elements in cmp 


B.2.12 get_sk 

Format: (sl,s2,s3^get_sk(k) 

Description. Extracts the switching matrices SHk\ S2(k) and S3(k ) for the k rt STP from the global workspace 
variables Slk, S2k. and S3k. See Appendix A. 

Inputs : k indicates the k! h STP of the current BTP. NOTE: k = 1 . 2 stppbtp 

Outputs: si . s2. s3 switching matrices for the k‘ h STP corresponding to SHk), S2(k) and S3(k) respectively 

Comments: This routine uses the global workspace variables generated by mropt_init 

B.2.13 make_noise 

Format : [r >= make^noise (rw, rv) 

Description: Creates a compound noise covariance matrix from rw and rv 

Inputs : rw, rv the process and sensor noise covariance matrices respectively 

Outputs : r compound noise matrix of the form 

rw 0 


r = 


0 rv 


I 


s 
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B.2.14 make_plt 


Format : 

If, g, w, h. v, nzcmp )= makeplt ( pltc, npll, ncmp, stp ) 


Description : 

Creates the compound plant matnces G„ W and V for the current plant penurbation. See 

Appendix A. 


/rtpute : 

pltc has the same form as pit in Table 3. 1 except it describes only a single plant penurbation 
nplt, ncmp, stp See Table 3.1. 


Outputs : 

f, g, w, h, v compound plant matrices F„ G„ Wand V for the cu rent plant penurbation respecuvely 

nzcmp the number states in the periodically time-varying representation of the mult.rate 
compensator 


B.2.15 mroptcheck 


Fonnat: 

mropt_check 


Description : 

mropt_check performs elemental error checking on the data in Table 3.2 and set the variable 
err--l if it finds an error in the data. mropt_check is called by mroptjnit. 


/rtpwte : 

None 


Outputs : 

None 


B.2.16 vnropt_globai 

l 

Format ; 

mropt_global 


Description : 

Called by mroptjnit to define 'he global workspace variables 


Inputs : 

None 

i 

Outputs : 

None 


B.2.17 mropt_fminu 

i 

* 

Format : 

[x, OPTIONS] = mropt_fminu ( FUN, x, OPTIONS. GRADFUN, Pi, P2, P3. P4, P 5 . 



P6. P7. P8. P9 P10) 


Description ; 

A modified version of FMINU from the Optimization.Toolbox. This version automatically reduces 
the search step sue when a destabilizing solution is encountered. 


/npatt : 

See FMINU 


Outputs : 

See FMINU 


Comments: 

Requires the Optimization_Toolbox 





APPENDIX C 


M-File mkm 


The following M-File creates the state space matrices for the Mass-Spring-Mass system of Section 4.0 

function (a,b, c,d] =mkm(r.l , m2 , k, b) 

% [a, b, c, d] =Kikm(ml ,m2 , k , b) 

if nargin-=4 %nominal plant 
ml = 1 ; 
m2 = . 1 ; 
k=l; 
b= 0.01; 
end; 

t=(l/ml)+(l/m2) ; 

a ~ [ 0 1 C 0 

0 0 k/ml b/ml 

0 0 0 1 

0 0 - k*t -b* t ] ; 

b= [0 0 

1/ml 0 
0 0 

-1 /ml 1 /m2 ] ; 

c=eye ( 4 ) ; 
d=zeros (4,2) ; 


Appendix d 

Script mropt_mkm 


The following script defines the workspace variables in Table 3.2 for the example problem of Section 4.0 
% Script: mropt_mkm.m 

% Example script for creating input data for mropt 

^ the following data required by nuropt_init 

% pit , nplt , wxx , rv, qla . q2a , ma 

% ucol , ncol , srow, crow 
% sz.su, sy.cmp, stp, Stppbtp, cmpfree 

% G. Mason, U.W. 1992 

% Continuous plant description 


Np=2 ; 
nplt=4 ,- 

Plt=[] ; 


% The number of plant conditions 

% The number of states in the plants 

% Clear the variable which holds tne plant data 


% Plant #1 

[a, b, c, d] =mkm ( 1 , 0.1, 1,0. 01) ; 
pit = [a b; c d] ; 

% Plant #2 

[a, b, c, d] =mkm (1,0.2, 1,0.01); 
pit = [pit [a b ; c a] ] ; 


% Pointer to rows and columns in the plant 

ucol=[l 2],- % control input 

ncol=2; % process noise input 

srow= [1 3],- % sensor output from the piant 

crow= ( 1 234]; % criterion output from the plant 

% Continuous process noise PSD 

% wxx is a compound matrix like pit, it contains wxx for plant #1 and #2 
wxx= [.1.1]; 


% Discrete sensor noise PSD 
% with same compound form as pit 
rv=eye ( 2 ) * 0 . 0 0 1 ; 
rv= [rv rv] ; 


^ ~ “ " Continuous Cost weighting matrices 

% Output criterion weights 

% Use le-8 instead of 0 so that qal & qa2 are positive definite 
% The synthesis algorithm would accept semi-definite qla v q2a 
'% but calc_byGcost will not 

qlal=diag( (5 . 5 le-8 le-8 2.2]); 
qla2=diag( [6.5 le-8 le-8 4.8]); 
qla= (qal/10 . 6 qa2/8.0]; 


% Control input weights 
q2a=diag( (le4 lei]); 
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q2a=[ra/10.6 ra/8.0]; 


% Cross weighting (ma) 

iS Cr ° SS " 9i 9hting it 


can be left undefined 


l 7“7" Compensator description 

cmp=? 0.5 P J OCesso J trains from successive loop closures design 

o 0 1 0 1 

0.08 0 - 0.2 0 
0 °- 6 o - 0 . 75 ] ; 


% Sampling schedule 

su= f n n . i ^ ^ . 


su=[0 0 ; 1 1 ; 0 0 ; 0 1 ; 0 0;0 1 ; 0 0;0 1] ; 

sy=[8 21; % "Z enSat ° r out P ut updating w/ delay 

sz= [ 3 2]* input sam Pling, multiplexed 

stp = .1- l Compensator state update 

Stppbtp = 8 - I Jlf Sh °^ test ^pling period 

% The nulnbec of stp' s in one BTP 

% Free compensator gains 
cmpfree=[l Oil 
0 111 
10 11 
0 1 1 1 ] • 


% clear temporary variables 
clear abed 


